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

2019-06-13 16:38:39 浏览数 (1)

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

代码语言：javascript复制
ls -lh *raw.vcf
代码语言：javascript复制
-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文件

代码语言：javascript复制
rm *_marked*.bam
代码语言：javascript复制
ls  *.bam  |xargs -i samtools index {} 
ls  *.bam  | while read id ;do (samtools flagstat $id > $(basename $id ".bam").stat);done
代码语言：javascript复制
cat SRR7696207.stat
代码语言：javascript复制
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

代码语言：javascript复制
conda install -c bioconda bedtools

制作exon.bed文件

代码语言：javascript复制
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"t0t "}'  > /mnt/f/kelly/bioTree/server/wesproject/align/hg38.exon.bed

查看

代码语言：javascript复制
cat hg38.exon.bed |head
代码语言：javascript复制
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文件夹

代码语言：javascript复制
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文件都移动到里面

代码语言：javascript复制
mkdir stats
mv *stat stats/
ls -lh stats/

显示如下

代码语言：javascript复制
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文件夹，把上面这种文件夹都移动到里面

代码语言：javascript复制
mkdir qualimap
mv *_stats qualimap
cd qualimap
ls -lh
代码语言：javascript复制
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

代码语言：javascript复制
multiqc ./

查看

multimap_multiqc

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

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

代码语言：javascript复制
multiqc ./

然后把得到的

代码语言：javascript复制
├── [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

marked

0 人点赞