用二代测序数据的reads组装一个基因序列

之前没有做过用二代测序数据的paired-end 数据组装一个基因。今天实验室有一个同学的在图位克隆的时候遇到了一个问题，发现有一个候选基因的可能性很大，从IGV浏览器中看到，这个基因在野生型材料和突变体材料之间有38个碱基的缺失。但是设计引物扩增的时候，无法扩增出来，于是想可能该38个碱基附近有很长的T-DNA序列的插入(这个突变体是EMS诱变的获得的， EMS诱变的那个材料也是一个突变体，我们怀疑是TDNA 插入的一个突变体），所以问我能否用这个区间的reads 组转出来。我觉得可以试一下。首先，我根据这个基因的位置信息，在比对的bam文件中，获取这些reads的ID，然后用这些ID在fastq 文件中获取这些reads的完整序列。因为我发现如果用Bam文件中的mapping上的reads 直接获取序列的话，有些reads的序列在比对的时候被截断了（bam文件中的soft clip和hard clip）。其次，获取了在区间的reads之后，用Abyss进行组装。不幸的是，这个地方可能复杂，没有组装成一个完整的contig。最后，我把参考基因组这段序列和这些contig进行了比较，发现都能比对到参考基因组上，但是有的contig的部分序列是比对不上的，这也就是无法组装成一个完整的contig的原因。问题是如果这个基因在野生型和突变体中，缺失38bp，怎么就扩增不出来呢？组装出的contig序列到底是否含有T-DNA 的序列吗？这个问题需要进一步的blast。未完，这个还要继续。

下边是分析一些代码：

from pysam import AlignmentFile
from pysam import  FastxFiledef write_read_id(bam, region, out_prefix):"""write reads ID in a given region in a sorted bam file"""out_fq = open(out_prefix+'.txt', 'w')sam = AlignmentFile(bam, 'rb')for read in sam.fetch(contig=region[0], start=region[1], end=region[2]):if read.is_paired:out_fq.write(read.query_name+"\n")out_fq.close()sam.close()def write_read_fastq(fastq_file, fq_id_file, out_prefix):""" look up ids from fa_id_file and write fastq format for ids"""out_fq = open(out_prefix+'.fastq', 'w')fq = FastxFile(filename=fastq_file)id_file = open(fq_id_file, 'r')ids = [item.strip() for item in id_file.readlines()]for read in fq:if read.name in ids:out_fq.write(str(read)+'\n')out_fq.close()fq.close()id_file.close()write_read_id('fpa-3.sorted.bam', region=('Chr5', 7741200, 7743400), out_prefix='Ath_ids')
write_read_fastq(fastq_file='fpa-3_1.fastq', fq_id_file='Ath_ids.txt', out_prefix='gene_fq_1')
write_read_fastq(fastq_file='fpa-3_2.fastq', fq_id_file='Ath_ids.txt', out_prefix='gene_fq_2')
#ubuntu
sudo apt-get install abyss
abyss-pe name=Ath_gene  k=25 B=1G in="gene_fq_1.fastq  gene_fq_2.fastq"

用二代测序数据的reads组装一个基因序列相关推荐

MPB：微生物所蔡磊组-基于二代测序的真菌基因组组装和注释
为进一步提高<微生物组实验手册>稿件质量,本项目新增大众评审环节.文章在通过同行评审后,采用公众号推送方式分享全文,任何人均可在线提交修改意见.公众号格式显示略有问题,建议电脑端点击文末阅 ...
关于文献中二代测序数据下载（NCBI）的问题
关于文献中二代测序数据下载(NCBI)的问题现在二代测序用于生物学研究非常广泛,大部分文章的序列会上传到Sequence Read Archive(SRA)上,这东西也属于NCBI数据库中的吧,我理 ...
三代测序数据超快组装软件--大牛Li heng 力作
三代测序数据超快组装软件--大牛Li heng 力作 (2017-06-19 16:53:46) 转载▼ 分类: 三代 1:软件链接:https://github.com/lh3/miniasm ...
用GATK进行二代测序数据 SNP Calling 流程：（二）bwa比对和HaplotypeCaller 变异检测
1. 创建基因组索引 bwa index genome.fa 2. 查看read group信息,按read group分组, 比对.合并,生成gvcf 由于数据太多,无法存储过多的中间文件,因此写了 ...
二代测序数据统计分析中为什么是负二项分布？
本文转载自"universebiologygirl",已获授权 1. 转录组数据统计推断的难题在RNA-seq中进行两组间的差异分析是最正常不过的了我们在其它实验中同样会遇到类 ...
用GATK进行二代测序数据 SNP Calling 流程：（四）变异过滤
GATK推荐的最好的过滤方式是用 VQSR功能,它通过机器学习算法来判断SNP的优劣,因此至少需要两个已存在的 SNP 数据集,一个是经过验证的高质量 SNP 数据集作为真集(如 HapMap),还需 ...
经典：基因组测序数据从头拼接或组装算法的原理
基因组测序数据的拼接/组装 (图片来源:google) 每一个物种的参考基因组序列(reference genome)的产生都要先通过测序的方法,获得基因组的测序读段(reads),然后再进行从头拼接 ...
iMeta | 青岛华大范广益组基于共标签测序数据的高质量宏基因组组装工具MetaTrass...
点击蓝字关注我们 MetaTrass:基于共标签测序数据的人类肠道微生物高质量宏基因组组装工具 https://doi.org/10.1002/imt2.46 RESEARCH ARTICLE ●2 ...
测序数据的处理方法及装置制造方法及图纸
测序数据的处理方法及装置制造方法及图纸技术编号:19389025阅读:109留言:0更新日期:2018-11-10 02:04 本发明专利技术公开了一种测序数据的处理方法及装置.其中,该方法包括:拆 ...
三代测序数据纠错的方法、装置和计算机可读存储介质与流程
三代测序数据纠错的方法.装置和计算机可读存储介质与流程文档序号:15616049发布日期:2018-10-09 21:24 导航: X技术> 最新专利>计算;推算;计数设备的制造及其应用 ...

用二代测序数据的reads组装一个基因序列

用二代测序数据的reads组装一个基因序列相关推荐

最新文章

热门文章