- android - 多次调用 OnPrimaryClipChangedListener
- android - 无法更新 RecyclerView 中的 TextView 字段
- android.database.CursorIndexOutOfBoundsException : Index 0 requested, 光标大小为 0
- android - 使用 AppCompat 时,我们是否需要明确指定其 UI 组件(Spinner、EditText)颜色
我有一个 txt 文件,它是一个转换后的 fasta 文件,其中只有一个我有兴趣分析的特定区域。看起来像这样
CTGGCCGCGCTGACTCCTCTCGCT
CTCGCAGCACTGACTCCTCTTGCG
CTAGCCGCTCTGACTCCGCTAGCG
CTCGCTGCCCTCACACCTCTTGCA
CTCGCAGCACTGACTCCTCTTGCG
CTCGCAGCACTAACACCCCTAGCT
CTCGCTGCTCTGACTCCTCTCGCC
CTGGCCGCGCTGACTCCTCTCGCT
我目前正在使用 Excel 对每个位置的核苷酸多样性进行一些计算。有些文件的读取次数约为 200,000 次,因此这使得 Excel 文件变得笨重。我想一定有一种更简单的方法可以使用 python 或 R 来做到这一点。
基本上,我想获取包含序列列表的 .txt 文件,并使用此方程 –p(log2(p)) 测量每个位置的核苷酸多样性。有谁知道除了 Excel 之外如何以其他方式完成此操作?
提前非常感谢您的帮助。
最佳答案
如果您可以从 fasta 文件工作,那可能会更好,因为有专门设计用于该格式的软件包。
在这里,我在 R 中给出了一个解决方案,使用包 seqinr
以及dplyr
(tidyverse
的一部分)用于操作数据。
如果这是您的 fasta 文件(基于您的序列):
>seq1
CTGGCCGCGCTGACTCCTCTCGCT
>seq2
CTCGCAGCACTGACTCCTCTTGCG
>seq3
CTAGCCGCTCTGACTCCGCTAGCG
>seq4
CTCGCTGCCCTCACACCTCTTGCA
>seq5
CTCGCAGCACTGACTCCTCTTGCG
>seq6
CTCGCAGCACTAACACCCCTAGCT
>seq7
CTCGCTGCTCTGACTCCTCTCGCC
>seq8
CTGGCCGCGCTGACTCCTCTCGCT
您可以使用 seqinr
包将其读入 R:
# 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
这是一个更直观的 View (使用 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 来可视化它,类似于我上面所做的,但也许使用 facets将每个文件的多样性分成不同的面板。
关于python - 从测序数据中解析信息,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/45281606/
我是一名优秀的程序员,十分优秀!