SILVA数据库全库下载序列的处理

处理用于QPCR引物设计的16S序列

最近在做肠道微生物的课题,搜索得知SILVA数据库是最近更新而且用的最多的,看网上的教程把其全库的序列下载了下来,没有比对的有200多兆,比对完的超过三个G,参考的那个微信公众号文章说只需要下载没有比对的,我还不信邪,把两个都下载下来了,一个解压后有3G,另一个有76g多,实在是难以处理,3g多的还勉强可以操作,于是就一小的文件做了筛选。

筛选用的是我刚入门的python,虽然水平挺菜,但是至少能用,水平也或许制约这我难以处理76g的文件。贴上我的一段筛选代码,及其简单,都没有什么复杂结构,水平啊!就是从中筛选出一个门或者属/种的16S序列。

fin = open('current_Bacteria_unaligned.fa')
fout = open('Actinobacteria.fasta', 'w')
phylum = ''
dict = set()
i = 1
flag = 0
for line in fin:
    if line.startswith('>'):
        flag = 0
        #phylum = line.strip().split(';')[-2]
        if 'Actinobacteria' in line:
        #if phylum == 'Akkermansia':
            flag = 1
            fout.write(line)
            #print(line, i)
            i += 1
        #break
        #dict.add(phylum)

    elif flag == 1:
        fout.write(line)
#print(dict)
print(i)
fout.close()

筛选完后发现太多了,特别是门,有几十万条,于是就再精选一下,把门中每个属只留一个,瞬间就变成千条左右了。

fin = open('Actinobacteria-filter.fasta')
fout = open('Actinobacteria-filter-1000.fasta', 'w')
phylum = ''
genusdict = {}
genus = ''
i = 0
flag = 0
for line in fin:
    if line.startswith('>'):
        flag = 0
        #if len(line.strip().split(';')) < 13:
        #   continue
        #else:
        genus = line.strip().split(';')[-2]
        #print(genus)
        if genus not in genusdict.keys():
            flag = 1
            genusdict[genus] = ''
            fout.write(line)
        else:
            continue
    elif flag == 1:
        fout.write(line)
    #break
print(len(genusdict))
fout.close()

好了,这样差不多就能愉快地比对找保守区域设计引物了。

菌门/属、种 总条数 筛选后(每属一条)
厚壁菌门 Firmicutes 775435 2642
变形菌门 Proteobacteria 936216 1121
拟杆菌门菌 Bacteroidetes 396435 316
放线菌门 Actinobacteria 362892 381
阿克曼氏菌属 Akkermansia 2777
普拉梭菌 Faecalibacterium prausnitzii** 340
解木聚糖拟杆菌 Bacteroides xylanisolvens** 7
瑞士乳杆菌 Lactobacillus helveticus** 1321
费氏丙酸杆菌 Propionibacterium freudenreichii** 18
厌氧球菌属 Anaerotruncus colihominis** 3
一致粪球菌 Coprococcus eutactus 2
Escherichia albertii** 78
普通拟杆菌 Bacteroides vulgatus** 35
脆弱拟杆菌 Bacteroides ovatus** 32
迟缓埃格特菌 Eggerthella lenta 28
罗斯氏菌属 Roseburia 5905

发表评论

电子邮件地址不会被公开。 必填项已用*标注