最近打算测一些干燥材料和标本材料的浅层基因组,有些材料降解严重,拼接的过程可能会比较困难。由于我们的目标主要是转录组的数据,所以试试看能不能在拼接之前将与转录组无关的序列直接在原始测序数据中删除。
minimap2和samtools安装都比较简单,这里就不行详细介绍了。
1.参考数据建库
首先我们选择一个转录组数据进行建库作为参考的比对序列
minimap2 -d 007.min 007_trinity.fasta
# 007.min为建立的索引库的名称, 007_trinity.fasta为准备的参考fasta序列文件
这样就建立了一个参考库
2.比对
minimap2 -ax map-ont 007.min 128_1.fq.gz -o 128_1.sam
# 007.min 是参考库, 128_1.fq.gz是需要筛选的基因组原始数据,128_1.sam 是输出的文件
minmap2 的常用参数
主要分成五大类,索引(Indexing),回帖(Mapping),比对(Alignment),输入/输出(Input/Output),预设值(Preset)。
-x :非常中要的一个选项,软件预测的一些值,针对不同的数据选择不同的值
map-pb/map-ont: pb或者ont数据与参考序列比对;
ava-pb/ava-ont: 寻找pd数据或者ont数据之间的overlap关系;
asm5/asm10/asm20: 拼接结果与参考序列进行比对,适合~0.1/1/5% 序列分歧度;
splice: 长reads的切割比对
sr: 短reads比对
-d :创建索引文件名
-a :指定输出格式为sa格式,默认为PAF
-Q :sam文件中不输出碱基质量
-R :reads Group信息,与bwa比对中的-R一致
-t:线程数,默认为3
-o: 输出文件
3使用samtools处理sam文件
sam文件(上一步得到的那种文件)是一种列表格式,用来记录reads比对到参考基因组上的信息,包括哪一条reads,比对到哪条基因组上的哪个位置,是一对一比对还是一对多比对,有无错配,错配是怎样的。因此就需要包含很多列的信息,下面具体展示sam格式。
第一列:是reads ID
第二列:是flag标记的总和
第三列:比对到参考序列上的染色体号。
第四列:为在参考序列上的位置
第五列:比对的质量值,MAPQ
第六列:代表比对结果的CIGAR字符串
第七列:mate比对到的染色体号,若是没有mate,则是*
第八列:比对到参考序列上的第一个碱基位置
第九列:Template的长度,
第十列:为read的序列
第十一列:为ASCII码格式的序列质量;
比对得到的sam或者bam不能直接用于下游的数据处理,需要进行很多处理。主要包括转换为bam,排序,合并不同lane的数据,对bam文件加头部信息等操作。而这些工作都可以使用samtools工具来进行操作。samtools顾名思义,是处理sam格式的工具合集。samtools主要包含以下几大功能:Indexing 建立索引,Editing 编辑文件,File operations,Statistics,统计相关功能;Viewing,查看,等
常用参数:
samtools选项参数:
-o:输出结果文件
-O:输出文件格式
-@:线程数目
--reference:参考序列
详细操作不讲了,直接提取需要的目标序列
samtools fasta -F 4 128_1.sam > 128_1.minimap2.fasta
# -F 4 参数表示提取出比对上的序列, -f 4 则是提取未比对上的序列。
最后得到一个fasta文件,使用这个fasta文件就可以进行下一步的序列拼接了。
之后使用Trinity进行拼接测试,但是总是运行一段时间就莫名其妙的断掉了,也没有任何报错信息,后来将输出的fasta文件改成fastq文件之后可以正常拼接,输出fastq的命令问
samtools fastq -F 4128_1.sam > 128_1.minimap2.fastq