使用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 secondary
372636 + 0 supplementary
0 + 0 duplicates
55374129 + 0 mapped (99.96% : N/A)
55026224 + 0 paired in sequencing
27513112 + 0 read1
27513112 + 0 read2
54512924 + 0 properly paired (99.07% : N/A)
54978908 + 0 with itself and mate mapped
22585 + 0 singletons (0.04% : N/A)
330146 + 0 with mate mapped to a different chr
252082 + 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;
do
sample=${id%%.*}
echo $sample
qualimap bamqc --java-mem-size=8G -gff hg38.exon.bed -bam $id &
done
align下新建stats文件夹，把stat文件都移动到里面
mkdir stats
mv *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 qualimap
mv *_stats qualimap
cd qualimap
ls -lh
total 0
drwxrwxrwx 0 root root 4.0K Jun 7 17:41 SRR7696207_bqsr_stats
drwxrwxrwx 0 root root 4.0K Jun 7 17:58 SRR8517853_bqsr_stats
drwxrwxrwx 0 root root 4.0K Jun 7 18:03 SRR8517854_bqsr_stats
drwxrwxrwx 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
