Если вы можете работать с файлом fasta, это может быть лучше, так как существуют пакеты, специально предназначенные для работы с этим форматом.
Здесь я даю решение на R, используя пакеты seqinr
, а также dplyr
(часть tidyverse
) для управления данными.
Если бы это был ваш файл fasta (на основе ваших последовательностей):
>seq1
CTGGCCGCGCTGACTCCTCTCGCT
>seq2
CTCGCAGCACTGACTCCTCTTGCG
>seq3
CTAGCCGCTCTGACTCCGCTAGCG
>seq4
CTCGCTGCCCTCACACCTCTTGCA
>seq5
CTCGCAGCACTGACTCCTCTTGCG
>seq6
CTCGCAGCACTAACACCCCTAGCT
>seq7
CTCGCTGCTCTGACTCCTCTCGCC
>seq8
CTGGCCGCGCTGACTCCTCTCGCT
Вы можете прочитать его в R, используя пакет seqinr
:
# Load the packages
library(tidyverse) # I use this package for manipulating data.frames later on
library(seqinr)
# Read the fasta file - use the path relevant for you
seqs <- read.fasta("~/path/to/your/file/example_fasta.fa")
Это возвращает объект list
, который содержит столько элементов, сколько последовательностей в вашем файле.
Для вашего конкретного вопроса - расчета показателей разнообразия для каждой позиции -
мы можем использовать две полезные функции из пакета seqinr
:
getFrag()
для подмножества последовательностей
count()
для вычисления частоты каждого нуклеотида
Например, если нам нужны частоты нуклеотидов для первой позиции наших последовательностей, мы могли бы сделать:
# Get position 1
pos1 <- getFrag(seqs, begin = 1, end = 1)
# Calculate frequency of each nucleotide
count(pos1, wordsize = 1, freq = TRUE)
a c g t
0 1 0 0
Показывая нам, что первая позиция содержит только «C».
Ниже приведен способ программного «перебора» всех позиций и выполнения вычислений, которые могут нас заинтересовать:
# Obtain fragment lenghts - assuming all sequences are the same length!
l <- length(seqs[[1]])
# Use the `lapply` function to estimate frequency for each position
p <- lapply(1:l, function(i, seqs){
# Obtain the nucleotide for the current position
pos_seq <- getFrag(seqs, i, i)
# Get the frequency of each nucleotide
pos_freq <- count(pos_seq, 1, freq = TRUE)
# Convert to data.frame, rename variables more sensibly
## and add information about the nucleotide position
pos_freq <- pos_freq %>%
as.data.frame() %>%
rename(nuc = Var1, freq = Freq) %>%
mutate(pos = i)
}, seqs = seqs)
# The output of the above is a list.
## We now bind all tables to a single data.frame
## Remove nucleotides with zero frequency
## And estimate entropy and expected heterozygosity for each position
diversity <- p %>%
bind_rows() %>%
filter(freq > 0) %>%
group_by(pos) %>%
summarise(shannon_entropy = -sum(freq * log2(freq)),
het = 1 - sum(freq^2),
n_nuc = n())
Результат этих вычислений теперь выглядит так:
head(diversity)
# A tibble: 6 x 4
pos shannon_entropy het n_nuc
<int> <dbl> <dbl> <int>
1 1 0.000000 0.00000 1
2 2 0.000000 0.00000 1
3 3 1.298795 0.53125 3
4 4 0.000000 0.00000 1
5 5 0.000000 0.00000 1
6 6 1.561278 0.65625 3
А вот более наглядный вид (используя ggplot2
, также часть пакета tidyverse
):
ggplot(diversity, aes(pos, shannon_entropy)) +
geom_line() +
geom_point(aes(colour = factor(n_nuc))) +
labs(x = "Position (bp)", y = "Shannon Entropy",
colour = "Number of\nnucleotides")
Обновление:
Чтобы применить это к нескольким файлам fasta, вот одна из возможностей (я не тестировал этот код, но что-то вроде этого должно работать):
# Find all the fasta files of interest
## use a pattern that matches the file extension of your files
fasta_files <- list.files("~/path/to/your/fasta/directory",
pattern = ".fa", full.names = TRUE)
# Use lapply to apply the code above to each file
my_diversities <- lapply(fasta_files, function(f){
# Read the fasta file
seqs <- read.fasta(f)
# Obtain fragment lenghts - assuming all sequences are the same length!
l <- length(seqs[[1]])
# .... ETC - Copy the code above until ....
diversity <- p %>%
bind_rows() %>%
filter(freq > 0) %>%
group_by(pos) %>%
summarise(shannon_entropy = -sum(freq * log2(freq)),
het = 1 - sum(freq^2),
n_nuc = n())
})
# The output is a list of tables.
## You can then bind them together,
## ensuring the name of the file is added as a new column "file_name"
names(my_diversities) <- basename(fasta_files) # name the list elements
my_diversities <- bind_rows(my_diversities, .id = "file_name") # bind tables
Это даст вам таблицу разнообразия для каждого файла. Затем вы можете использовать ggplot2
для визуализации, аналогично тому, что я сделал выше, но, возможно, используя фасеты разделить разнообразие из каждого файла на разные панели.
person
hugot
schedule
25.07.2017
dput()
ваши данные с ожидаемыми результатами. - person Christoph   schedule 24.07.2017