8 比对及找变异步骤的质控

使用qualimap对wes的比对bam人家总结测序深度和覆盖度

ls -lh *raw.vcf

-rwxrwxrwx 1 root root 184M Jun 7 10:58 SRR7696207_raw.vcf-rwxrwxrwx 1 root root 61M Jun 7 09:39 SRR8517853_raw.vcf-rwxrwxrwx 1 root root 87M Jun 7 03:04 SRR8517854_raw.vcf-rwxrwxrwx 1 root root 331M Jun 7 02:21 SRR8517856_raw.vcf

1 比对的各个阶段的bam进行质控

可以把中间生成的.bam文件删除，就是带marked的bam文件

rm *_marked*.bam

ls *.bam |xargs -i samtools index {} ls *.bam | while read id ;do (samtools flagstat $id > $(basename $id ".bam").stat);done

cat SRR7696207.stat

55398860 + 0 in total (QC-passed reads + QC-failed reads)0 + 0 secondary372636 + 0 supplementary0 + 0 duplicates55374129 + 0 mapped (99.96% : N/A)55026224 + 0 paired in sequencing27513112 + 0 read127513112 + 0 read254512924 + 0 properly paired (99.07% : N/A)54978908 + 0 with itself and mate mapped22585 + 0 singletons (0.04% : N/A)330146 + 0 with mate mapped to a different chr252082 + 0 with mate mapped to a different chr (mapQ>=5)

安装bedtools

conda install -c bioconda bedtools

制作exon.bed文件

cat /mnt/f/kelly/bioTree/server/wesproject/hg38/annotation/CCDS.20160908.txt |perl -alne '{/\[(.*?)\]/;next unless $1;$gene=$F[2];$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$_\t$gene" foreach split/,/,$exons;}'|sort -u |bedtools sort -i |awk '{print "chr"$0"\t0\t+"}' > /mnt/f/kelly/bioTree/server/wesproject/align/hg38.exon.bed

查看

cat hg38.exon.bed |head

chr1 69090 70007 OR4F5 0 +chr1 450739 451677 OR4F29 0 +chr1 685715 686653 OR4F16 0 +chr1 925941 926012 SAMD11 0 +chr1 930154 930335 SAMD11 0 +chr1 931038 931088 SAMD11 0 +chr1 935771 935895 SAMD11 0 +chr1 939039 939128 SAMD11 0 +chr1 939274 939459 SAMD11 0 +chr1 941143 941305 SAMD11 0 +

qualimap进行质控

align文件夹

ls *_bqsr.bam | while read id;dosample=${id%%.*}echo $samplequalimap bamqc --java-mem-size=8G -gff hg38.exon.bed -bam $id & done

align下新建stats文件夹，把stat文件都移动到里面

mkdir statsmv *stat stats/ls -lh stats/

显示如下

total 0-rwxrwxrwx 1 root root 453 Jun 7 16:31 SRR7696207_bqsr.stat-rwxrwxrwx 1 root root 447 Jun 7 16:29 SRR7696207.stat-rwxrwxrwx 1 root root 444 Jun 7 16:34 SRR8517853_bqsr.stat-rwxrwxrwx 1 root root 444 Jun 7 16:33 SRR8517853.stat-rwxrwxrwx 1 root root 447 Jun 7 16:37 SRR8517854_bqsr.stat-rwxrwxrwx 1 root root 447 Jun 7 16:35 SRR8517854.stat-rwxrwxrwx 1 root root 452 Jun 7 16:43 SRR8517856_bqsr.stat-rwxrwxrwx 1 root root 452 Jun 7 16:40 SRR8517856.stat

完成后会生成SRR8517856_bqsr_stats类似的文件夹

现在建立一个qualimap文件夹，把上面这种文件夹都移动到里面

mkdir qualimapmv *_stats qualimapcd qualimapls -lh

total 0drwxrwxrwx 0 root root 4.0K Jun 7 17:41 SRR7696207_bqsr_statsdrwxrwxrwx 0 root root 4.0K Jun 7 17:58 SRR8517853_bqsr_statsdrwxrwxrwx 0 root root 4.0K Jun 7 18:03 SRR8517854_bqsr_statsdrwxrwxrwx 0 root root 4.0K Jun 7 17:41 SRR8517856_bqsr_stats

然后做multiqc

multiqc ./

查看

multimap_multiqc

coverage不够，不知是我操作哪步有问题还是？

然后在stats文件夹下执行multiqc命令

multiqc ./

然后把得到的

├── [4.0K] multiqc_data│ ├── [ 261] multiqc_general_stats.txt│ ├── [7.3K] multiqc.log│ ├── [2.2K] multiqc_samtools_flagstat.txt│ └── [ 882] multiqc_sources.txt├── [1.0M] multiqc_report.html

下载到本地电脑查看。

stats_multiqc