计算paired-end测序的insert size

给出的条件为:
参考序列文件: 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数会偏少。

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据