我对一个文件很感兴趣。看起来像这样
CTGGCCGCTGCTCCTCTCGCT公司
CTCGCAGCATCTCCTCTTGCG公司
ctagcgcctctgactccgctagcg
CTCGCTGCCCTCTCACCTCTTGCA
CTCGCAGCATCTCCTCTTGCG公司
CTCGCAGCACTAACCCTAGTAGCT
CTCGCTGCTCTGCTCCTCGCC
CTGGCCGCTGCTCCTCTCGCT公司
我目前正在使用excel对每个位置的核苷酸多样性进行一些计算。有些文件有大约200000次读取,所以这使得excel文件难以处理。我想使用python或R一定有更简单的方法来实现这一点
基本上,我想用带有序列列表的.txt文件,用这个等式-p(log2(p))测量每个位置的核苷酸多样性。有人知道除了excel之外,还有什么方法可以做到这一点吗?在
非常感谢您的帮助。在
您可以打开并读取文件:
然后用你
plist
来计算你的熵如果你能从档案里找到更好的方法 专门设计用于该格式的包。在
在这里,我在R中使用包
seqinr
和dplyr
(是tidyverse
)的一部分,用于操作数据。在如果这是你的fasta文件(根据你的顺序):
您可以使用
^{pr2}$seqinr
包将其读入R:这将返回一个
list
对象,该对象包含的元素数量与现有的一样多 文件中的序列。在对于您的特定问题-计算每个职位的多样性指标-
我们可以使用
seqinr
包中的两个有用函数:getFrag()
对序列进行子集化count()
计算每个核苷酸的频率例如,如果我们想知道 我们的序列,我们可以:
告诉我们第一个位置只包含一个“C”。在
下面是一种通过编程“循环”所有位置并执行 我们可能感兴趣的计算:
这些计算的输出如下所示:
下面是一个更直观的视图(使用
ggplot2
,也是tidyverse
包的一部分):更新:
要将其应用于多个fasta文件,有一种可能 (我没有测试此代码,但类似这样的代码应该可以工作):
这将为您提供每个文件的多样性表。然后,您可以使用
ggplot2
来可视化它,类似于我上面所做的,但是也许可以使用facets将每个文件的多样性分离到不同的面板中。在相关问题 更多 >
编程相关推荐