碱基编辑二代测序数据之筛选序列
生信人R语言学习必备
立刻拥有一个Rstudio账号
开启升级模式吧
（56线程，256G内存，个人存储1T）
一般来说二代测序数据的数据量是非常庞大的，其中也包含了一些我们不需要的数据，经过fastp软件筛选之后，还存在大量不符合我们期望的数据，那如何进一步筛选我们想要的数据呢，今天就让小果带大家看看数据是怎么被进一步筛选的吧。
下面小果将用碱基编辑的二代测序数据做一个详细的步骤介绍。首先我们要拥有碱基编辑前的原序列信息，这样才能对我们编辑后的信息做对比，原信息如图所示。
import pandas as pdflank = 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 = 0o = 0k = 0for 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].endgN20_right_start = gN20_right_index[0].startgn20 = line[gN20_left_end:gN20_right_start]#确定GN20if gn20 in GN20_list:o+=1N20_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 positionN20_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-29n20 = line[N20_left_end+4:N20_right_start-3] #确定N20if similarity(gn20,n20)>0.7:#相似度N20.write(gn20+','+n20+','+line)k+=1if i%500000==0:#程序进程print(sample_barcode,i,o,k)i+=1end = time.time()print(end-start)#记录时间B4.close()N20.close()
当脚本运行结束，我们就得到了符合条件的所有序列。
好了，今天小果的分享就到这里啦，小伙伴们有什么问题欢迎来和小果讨论分享呀。
如果小伙伴有其他数据分析需求，可以尝试使用本公司新开发的生信分析小工具云平台，零代码完成分析，非常方便奥，云平台网址为：(http://www.biocloudservice.com/home.html)，其中也包括了通路表达分析(http://www.biocloudservice.com/313/313.php)，单细胞的基因共表达分析(http://www.biocloudservice.com/906/906.php)等各种小工具哦~，有兴趣的小伙伴可以登录网站进行了解。
点击“阅读原文”立刻拥有
↓↓↓