2024-06-16

二代测序数据根据barcode分组

原创小果生信果 2023-08-16 19:00:43

生信人R语言学习必备

立刻拥有一个Rstudio账号

开启升级模式吧

（56线程，256G内存，个人存储1T）

还在为如何入门生信而彷徨吗？还在为Linux系统而感到无助吗？快关注小果的微信公众号，小果手把手带你入门生信，没有门槛、没有难度，跟着效果走，啥问题都没有。

当我们拿到的二代测序数据的数据量非常庞大时，我们怎么对其进行分组呢，今天小果来带你探索其中的奥妙。

二代测序数据经过fastp软件筛选之后我们就得到了成千上万条如图所示的序列，每组序列我们都在它的前端添加相同的barcode如图中红色标记所示，这样我们就可以用遍历的方式将相同组的序列放在一起。

因为双端测序存在反向互补的问题，所以首先我们先定义一个反向互补的函数。

#互补函数def complement(sequence):    sequence = sequence.upper()    sequence = sequence.replace('A', 't')    sequence = sequence.replace('T', 'a')    sequence = sequence.replace('C', 'g')    sequence = sequence.replace('G', 'c')    return sequence.upper()#反向函数def reverse(sequence):    sequence = sequence.upper()    return sequence[::-1]#反向互补函数def complement_reverse(sequence):return complement(reverse(sequence))

函数编写完成后，我们就可以将barcode放入函数中，就可以输出反向互补的序列啦。

print(complement_reverse('CCATCAATGCC') ,complement_reverse('CGATGGCGATA'),complement_reverse('ACACAAGCACC'),complement_reverse('CCGTTTCGACG'),complement_reverse('GGAAGTAGACC'))

根据barcode分组

rf = open("./输出文件.fq") #改成自己二代测序的fq文件A1= open("./1_A1.txt",'w')#创建txt文件写入不同组数据A2= open("./1_A2.txt",'w')A3= open("./1_A3.txt",'w')A4= open("./1_A4.txt",'w')A5= open("./1_A5.txt",'w')i=0for line in rf:  #遍历二代测序数据    line = line.strip()    if ('CCATCAATGCC' in line):#若序列存在CCATCAATGCC，则写入A1文件中        A1.write(line)        A1.write('n')    elif ('GGCATTGATGG'in line):        A1.write(complement_reverse(line))        A1.write('n')
    elif ('CGATGGCGATA' in  line):        A2.write(line)        A2.write('n')    elif ('TATCGCCATCG' in line):        A2.write(complement_reverse(line))        A2.write('n')
    elif ('ACACAAGCACC' in  line):        A3.write(line)        A3.write('n')    elif ('GGTGCTTGTGT' in  line):        A3.write(complement_reverse(line))        A3.write('n')
    elif ('CCGTTTCGACG' in  line):        A4.write(line)        A4.write('n')    elif ('CGTCGAAACGG' in  line):        A4.write(complement_reverse(line))        A4.write('n')            elif ('GGAAGTAGACC' in  line):        A5.write(line)        A5.write('n')    elif ('GGTCTACTTCC' in  line):        A5.write(complement_reverse(line))        A5.write('n')
    i+=1#     if i >10000:#         break    if (i%10000000 ==0): #查看代码进度        print(i)rf.close() #关闭文件