给出的条件为:
参考序列文件: ref.fasta
Illumnina paired-end文件: reads1.fq, reads2.fq
1. 结合 Bowtie2,SAMtools 和 Picard的使用来计算insert size
1.1 将paired-end数据比对到参考序列上
$ bowtie2-build ref.fasta ref $ bowtie2 --rg-id librarylength_insert --rg "PL:ILLUMINA" \ --rg "SM:Strain" -x ref -p 8 \ -1 reads1.fq -2 reads2.fq -S align.sam
1.2 将sam格式比对结果转换为bam格式并排序
使用samtools将sam文件转换成bam文件;使用picard将bam文件进行排序.
$ samtools view -bS align.sam > align.bam $ jar Xmx2g -jar $picardHome/SortSam.jar \ I=align.bam O=align.sort.bam SO=coordinate
1.3 计算insertsize
输入排序后的bam文件,使用picard生成insertsize的txt的结果文件和PDF格式的图形结果。
$ jar Xmx2g -jar $picardHome/CollectInsertSizeMetrics.jar \ I=align.sort.bam R=ref.fasta \ O=insertsize.txt H=insertsize.pdf
2. 使用Qualimap来计算insertsize
2.1 Qualimap的介绍
Qualimap 用于NGS比对数据的质量控制。使用JAVA和R编写的程序,有GUI和command-line两种运行方式。
文献:Fernando García-Alcalde, Konstantin Okonechnikov, José Carbonell, Luis M. Cruz, Stefan Götz, Sonia Tarazona, Joaquín Dopazo, Thomas F. Meyer, and Ana Conesa “Qualimap: evaluating next-generation sequencing alignment data.” Bioinformatics 28, no. 20 (2012): 2678-2679.
doi: 10.1093/bioinformatics/bts503
2.2 Qualimap的安装
Qualimap的正常使用需要 JAVA runtime vesion 6 或以上版本; R 2.14 或以上版本;需要一些 R 包。该软件默认使用最高线程数运行,速度不错。
$ wget http://qualimap.bioinfo.cipf.es/release/qualimap_v0.7.1.zip $ unzip qualimap_v0.7.1.zip $ cd qualimap_v0.7.1 $ sudo su $ /usr/local/bin/Rscript scripts/installDependencies.r $ R > source("http://bioconductor.org/biocLite.R") > biocLite("gplots","Repitools")
3. 注意事项
对于基因组的paired-end测序,将reads比对到参考基因组上;但是对于转录组测序的paired-end测序,则最好将reads比对到转录组序列上。而比对到基因组序列上reads数会偏少。