专栏首页科技记者SILVA、GREENGENES、RDP三大数据库的序列探索统计

SILVA、GREENGENES、RDP三大数据库的序列探索统计

最近对16s的三大数据库的序列的具体序列情况挺好奇的,决定统计一下各个序列的长度分布情况,以及这些序列具体分布在哪几个V区,有助于我解决后面16So数据的问题。还是用上我三脚猫的功夫,开始今天 的探索,没有人探索的事情,还是挺开心的。

1.统计序列长度分布情况

01
#获得长度列表文件
02
length_list = []
03
with open('current_Bacteria_unaligned.fa') as f:
04
    flag = 0
05
    length = 0
06
    for line in f:
07
        if line.startswith('>'):
08
            flag = 1
09
            if length != 0:
10
                #print(length)
11
                length_list.append(length)
12
                #break
13
                length = 0
14
            else:
15
                length = 0
16
        elif flag == 1:
17
            length +=len(line)
18

19
fout =  open('length_table.txt', 'w')
20
for a in length_list:
21
    fout.write(str(a) + '\n')
22
fout.close()
01
#统计长度区间分布
02
length_150 = 0
03
length_300 = 0
04
length_600 = 0
05
length_1300 = 0
06
length_1300_ = 0
07
with open('length_table.txt') as f:
08
    for line in f:
09
        if 0 < int(line.strip()) <= 150:
10
            length_150 += 1
11
        elif 150 < int(line.strip()) <= 300:
12
            length_300 += 1
13
        elif 300 < int(line.strip()) <= 600:
14
            length_600 += 1
15
        elif 600 < int(line.strip()) <= 1300:
16
            length_1300 += 1
17
        elif int(line.strip()) > 1300:
18
            length_1300_ += 1
19
print('150以内:', length_150,'\n',
20
 '150-300:', length_300, '\n',
21
  '300-600:', length_600, '\n',
22
  '600-1300:', length_1300,'\n',
23
  '1300-:', length_1300_)
1
#最后结果
2
150以内: 0
3
 150-300: 0
4
 300-600: 617,554
5
 600-1300: 1,063,377
6
 1300-: 1,515,109

我也可视化一下:

更正一下,这里用的是RDP数据库,之前由于Silva数据库用的是不兼容的14级分类系统而没采用。

从图中可以看出,大部分序列还是集中在600以上。

接着是greengenes数据库,这个数据库虽然序列较少,但是长度大部分集中在1300+,质量较高,就是好久没更新过了。 150以内: 0 150-300: 0 300-600: 0 600-1300: 895 1300-: 1262090

2.统计V区分布情况

从一个公众号得到的一张分布图是这样的,

我想确定的是序列都包含在哪个或者哪两个区。花了好大有力气,才把代码理好,效率应该还是比较低的,但是,对于我来说

001
import re
002

003
dic ={}
004
v1f = re.compile('GGATCCAGACTTTGATYMTGGCTCAG', re.I)
005
v3f = re.compile('CCTA[CT]GGG[AG][GTC]GCA[CG]CAG', re.I)
006
v4f = re.compile('GTG[CT]CAGC[AC]GCCGCGGTAA', re.I)
007
v4r = re.compile('ATTAGA[AT]ACCC[CTG][ATGC]GTAGTCC', re.I)
008
v6f = re.compile('AACGCGAAGAACCTTAC', re.I)
009
v8f = re.compile('CGTCATCC[AC]CACCTTCCTC', re.I)
010
vr = re.compile('AAGTCGTAACAAGGTA[AG]CCGTA', re.I)
011
#v4r = re.compile('GGACTAC[ATGC][ACG]GGGT[AT]TCTAAT', re.I)
012
vf = re.compile('AG[AG]GTT[CT]GAT[CT][AC]TGGCTCAG', re.I)
013
#vr = re.compile('GACGGGCGGTG[AT]GT[AG]CA', re.I)
014
#seq_list = []
015

016

017
def decide_which_zone(f1, f3, f4, f4r, f6, f8, f, fr):
018
    type = [f1, f3, f4, f4r, f6, f8, f, fr]
019
    type_str = ['f1', 'f3', 'f4', 'f4r', 'f6', 'f8', 'f', 'fr']
020
    type_name = []
021
    for i in range(len(type)):
022
        #print(type[i])
023
        if type[i] != None :
024
            type_name.append(type_str[i])
025
        else:
026
            continue
027
        if i + 1 < len(type):
028
            if type[i + 1] != None and type_str[i + 1] not in type_str:
029
                type_name.append(type_str[i + 1])
030
        else:
031
            continue
032
        if i + 2 < len(type):
033
            if type[i + 2] != None and type_str[i + 2] not in type_str:
034
                type_name.append(type_str[i + 2])
035
        else:
036
            continue
037
        if  i + 3 < len(type):
038
            if type[i + 3] != None and type_str[i + 3] not in type_str:
039
                type_name.append(type_str[i + 3])
040
        else:
041
            continue
042
        if  i + 4 < len(type):
043
            if type[i + 4] != None and type_str[i + 4] not in type_str:
044
                type_name.append(type_str[i + 4])
045
        else:
046
            continue
047
        if  i + 5 < len(type):
048
            if type[i + 5] != None and type_str[i +5 ] not in type_str:
049
                type_name.append(type_str[i + 5])
050
        else:
051
            continue
052
        if  i + 6 < len(type):
053
            if type[i + 6] != None and type_str[i + 6] not in type_str:
054
                type_name.append(type_str[i + 6])
055
        else:
056
            continue
057
        if  i + 7 < len(type):
058
            if type[i + 7] != None and type_str[i + 7] not in type_str:
059
                type_name.append(type_str[i + 7])
060
        else:
061
            continue
062
    return type_name
063

064
with open('current_Bacteria_unaligned.fa') as f:
065
    li = ['f1', 'f3', 'f4', 'f4r', 'f6', 'f8', 'f', 'fr']
066
    for i in range(len(li)):
067
        for j in range(len(li)):
068
            if i <= j:
069
                dic[li[i] + li[j]] = 0
070

071

072
    #print(dic)
073
    flag = 0
074
    seq = ''
075
    #i = 0
076
    for line in f:
077
        if line.startswith('>'):# and i <=2:
078
            flag = 1
079
            if seq != '':
080
                #print(seq)
081
                f1 = v1f.search(seq)
082
                f4r = v4r.search(seq)
083
                f3 = v3f.search(seq)
084
                f4 = v4f.search(seq)
085
                f6 = v6f.search(seq)
086
                f8 = v8f.search(seq)
087
                f = vf.search(seq)
088
                fr = vr.search(seq)
089
                #seq_list.append(length)
090
                type_name = []
091
                type_name = decide_which_zone(f1, f3, f4, f4r, f6, f8, f, fr)
092
                print(type_name)
093
                if type_name != []:
094
                    type_key = type_name[0] + type_name[-1]
095
                    dic[type_key] += 1
096
                #i += 1
097
                #break
098
                seq = ''
099
                flag = 0
100
            else:
101
                seq = ''
102
        elif flag == 1:
103
            seq += line.strip()
104
print(dic)
105
#运行结果
106
{'f1f1': 0, 'f1f3': 0, 'f1f4': 1, 'f1f4r': 0, 'f1f6': 0, 'f1f8': 0, 'f1f': 0, 'f1fr': 0, 'f3f3': 120867, 'f3f4': 143686, 'f3f4r': 356027, 'f3f6': 400108, 'f3f8': 4, 'f3f': 171637, 'f3fr': 53716, 'f4f4': 17211, 'f4f4r': 71719, 'f4f6': 40223, 'f4f8': 4, 'f4f': 13185, 'f4fr': 3646, 'f4rf4r': 39537, 'f4rf6': 45160, 'f4rf8': 0, 'f4rf': 977, 'f4rfr': 4263, 'f6f6': 43983, 'f6f8': 0, 'f6f': 64, 'f6fr': 2300, 'f8f8': 0, 'f8f': 0, 'f8fr': 0, 'ff': 1670, 'ffr': 26, 'frfr': 3678}

上图依然是RDP数据库,前面标错了。应该是某个地方出了问题,序列不全,少了估计有一半。不过趋势还是比较明显的,V3-V4是最多的,这也是前些年主流的研究的测序方式。

本文分享自微信公众号 - 科技记者(kejijizhe),作者:zd200572

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2018-08-21

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • SNP2HLA之参考数据集合并提高分型准确性

    这几天随便搜索snp2hla软件的参考数据集的时候发现一个韩国科学家写了一个数据集合并脚本,在使用韩国人样本测试时准确性较分别只用两个未合并的数据集准确性有所提...

    用户1075469
  • 肠道菌群学习笔记

    是指发酵糖类主要产物为乳酸的一类无芽孢、革兰氏染色阳性细菌的总称。它不是一个分类学上的概念,至少包含5种23属,其具有种和菌株的决定性,21种可用于食品和保健品...

    用户1075469
  • ubiome类似数据dada2处理探索2

    首先把usearch申请下载到工作目录,然后docker挂载到home,当然如果linux就直接省了这一步了,可以下载(安装)好直接使用。

    用户1075469
  • Go语言中反射的正确使用

    介绍 反射是元数据编程的一种形式,指的是程序获得本身结构的一种能力。不同语言的反射模型实现不一样,本文中的反射,仅仅指的是Go语言中的反射模型。 反射有两个问题...

    李海彬
  • Python面试题之Python中type和object的关系

    下面是jeff kit的回答: 给别人讲解过很多次,但写成文字是第一次。试一试吧,自己主要也是看了这篇文章(Python Types and Objects...

    Jetpropelledsnake21
  • Go语言中反射的正确使用

    介绍 反射是元数据编程的一种形式,指的是程序获得本身结构的一种能力。不同语言的反射模型实现不一样,本文中的反射,仅仅指的是Go语言中的反射模型。 反射有两个问题...

    李海彬
  • Go语言中反射的正确使用

    介绍 反射是元数据编程的一种形式,指的是程序获得本身结构的一种能力。不同语言的反射模型实现不一样,本文中的反射,仅仅指的是Go语言中的反射模型。 反射有两个问题...

    李海彬
  • oracle中的数组(第一篇)(r4笔记第9天)

    数组在各种编程语言中都是很重要的数据结构实现,在oracle中也有自己的一席之地。自己简单做了几个实验,发现很多东西还是眼高手低,真实去做的时候,里面还是有不少...

    jeanron100
  • dotnet检测类型是否为泛型

    private static string GetTableName(Type type) { //检测类型是否为泛型...

    hbbliyong
  • Json Schema 快速入门

    Json schema 本身遵循Json规范,本身就是一个Json字符串,先来看一个例子

    软测小生

扫码关注云+社区

领取腾讯云代金券