二代测序数据筛选序列
一般来说二代测序数据的数据量是非常庞大的，其中也包含了一些我们不需要的数据，经过fastp软件筛选之后，还存在大量不符合我们期望的数据，那如何进一步筛选我们想要的数据呢，今天就让小云带大家看看数据是怎么被进一步筛选的吧。
下面小云将用碱基编辑的二代测序数据做一个详细的步骤介绍。首先我们要拥有碱基编辑前的原序列信息，这样才能对我们编辑后的信息做对比，原信息如图所示。
import pandas as pd
flank = pd.read_csv("./TCDG_TDGNG_CDGNGG-第三版.csv", index_col=None, header=0,sep = ',',encoding="utf-8")
flank
拥有序列原信息后，就能从编辑后的序列中挑选我们想要的序列了，我们拿五组已经经过barcode分组的序列信息，根据barcode分组的详细步骤已经在上一篇文章中讲述，还没有掌握的小伙伴记得去看哦，我们将这五组的数据进行筛选。
筛选的原理就是我们利用GN20与N20的上下游序列不变，来根据上下游序列定位出GN20与N20，将定位出来的信息与原表格信息进行对比，我们规定如果对比相似度在0.7以上则留下该序列。
from Bio.Seq import Seq#导入库
from fuzzysearch import find_near_matches#序列对比
from difflib import SequenceMatcher#导入库
def similarity(a, b):#相似度对比
return SequenceMatcher(None, a, b).ratio()
import time#计时
GN20_list = list(flank['GN20'])#提取表格信息
for sample_barcode in ['A1','A2','A3','A4','A5']:#遍历每个txt文件
start = time.time()
B4 = open("./1_%s.txt"%sample_barcode)
N20 = open('./1_%s_N20.txt'%sample_barcode,'w')#创建新的txt，写入结果
i = 0
o = 0
k = 0
for line in B4:#遍历每条序列
gN20_right = 'GTTTTAGAGCTAGAAATA'
N20_right = 'AAAGAATTCTCGACCT'
gN20_left = 'GTGGAAAGGACGAAACACC'
gN20_left_index = find_near_matches(gN20_left, line, max_l_dist=2)#序列对比
if gN20_left_index:
gN20_right_index = find_near_matches(gN20_right, line, max_l_dist=2)
if gN20_right_index:
gN20_left_end = gN20_left_index[0].end
gN20_right_start = gN20_right_index[0].start
gn20 = line[gN20_left_end:gN20_right_start]#确定GN20
if gn20 in GN20_list:
o+=1
N20_right_index = find_near_matches(N20_right, line, max_l_dist=2)
if N20_right_index:
N20_right_start = N20_right_index[0].start
#ACGT relative position
N20_left_index = line[N20_right_start-29:N20_right_start-17].rfind('ACGT')
if N20_left_index>=0:
N20_left_end = N20_left_index + N20_right_start-29
n20 = line[N20_left_end+4:N20_right_start-3] #确定N20
if similarity(gn20,n20)>0.7:#相似度
N20.write(gn20+','+n20+','+line)
k+=1
if i%500000==0:#程序进程
print(sample_barcode,i,o,k)
i+=1
end = time.time()
print(end-start)#记录时间
B4.close()
N20.close()
当脚本运行结束，我们就得到了符合条件的所有序列。
好了，今天小云的分享就到这里啦，小伙伴们有什么问题欢迎来和小云讨论分享呀。
