专栏首页Y大宽8 比对及找变异步骤的质控

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 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

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 2 服务器基本情况

    Y大宽
  • Linux下查看进程的启动和运行时间

    总体来说, ps主要是查看进程的,尤其你关心的进程 top主要看cpu,内存使用情况,及占用资源最多的进程由高到低排序,关注点在于资源占用情况

    Y大宽
  • 7 对比对和变异结果用IGV进行可视化

    把上述文件下载到本地IGV查看 注意,igv同时需要.bam和相应的.bai文件,所以需要把整个文件夹cp。

    Y大宽
  • 如何使用CP / SCP / RSYNC在Linux中排除特定目录?

    对于任何系统管理员或一般Linux操作系统用户而言,在服务器之间执行文件复制操作都是一项常见任务。在将文件从一个系统复制到另一个系统时,由于某些特定原因,我们可...

    用户6543014
  • 15 道二叉树手写算法题(二)

    在上一期讲到,树和链表的手写算法题在面试中出现的频率最高。也正是因为这样,如果你马上就要参加面试,但之前没有刷多少算法题,那么很建议你先看看树和链表相关的题目。...

    乔戈里
  • admin3

    #################################################### 真机上实现别名的定义,修改配置文件

    py3study
  • CentOS7.4中Docker以rw方式挂载volume报Permission denied的解决思路

    版权声明:本文为耕耘实录原创文章,各大自媒体平台同步更新。欢迎转载,转载请注明出处,谢谢

    耕耘实录
  • Minimum Distance Between BST Nodes

    用户1147447
  • 226 Invert Binary Tree

    /** * Definition for a binary tree node. * function TreeNode(val) { * thi...

    用户1624346
  • 剑指Offer - 面试题54. 二叉搜索树的第k大节点(二叉树循环遍历)

    来源:力扣(LeetCode) 链接:https://leetcode-cn.com/problems/er-cha-sou-suo-shu-de-di-kd...

    Michael阿明

扫码关注云+社区

领取腾讯云代金券