5 一步法找基因变异流程

1 先查看sam文件

随机选择3个

$ samtools mpileup SRR8517854.bam |head -95|tail -3
[mpileup] 1 samples in 1 input files
chr1    10105   N       8       AAAAcAAA        kuuu>6uK
chr1    10106   N       10      CCCCcCCCC^$c    uuuukAuuAK
chr1    10107   N       10      CCCCcCCCCc      upukaFfuKF

再看另外一个:

$ samtools mpileup SRR8517854.bam |head -168944|tail -3
[mpileup] 1 samples in 1 input files
chr1    909515  N       1       c       K
chr1    909516  N       1       c       K
chr1    909517  N       1       g       K

关于samtools mpileup的具体用法和结果解释请看

2 最简单的找变异流程

align文件夹

ref=/mnt/f/kelly/bioTree/server/wesproject/hg38/Homo_sapiens_assembly38.fasta
time samtools mpileup -ugf $ref *.bam|bcftools call -vmO z -o wxs_out.vcf.gz
less -S wxs_out.vcf.gz

3载入IGV查看

3.1先构建索引

批量构建

ls *.bam|xargs -i samtools index {}

如果是服务器需要下载到本地查看

4 去除PCR重复

需要samtools的4个命令

align文件夹

samtools markdup -r SRR7696207.bam SRR7696207.rm.bam
[markdup] error: no ms score tag. Please run samtools fixmate on file first.
[markdup] error: no ms score tag. Please run samtools fixmate on file first.

提示先fixmate

$ samtools fixmate SRR7696207.bam SRR7696207.fixmate.bam
[bam_mating_core] ERROR: Coordinate sorted, require grouped/sorted by queryname.

提示先要sort,并且要以queryname进行sort

$ samtools sort -n -o SRR7696207.sort.bam SRR7696207.bam
$ samtools sort -n -o SRR7696207.namesort.bam SRR7696207.bam
[bam_sort_core] merging from 25 files and 1 in-memory blocks...

接下来就可以逆回去了,也就是是正确的顺序(本篇最下方有说明)

$ samtools sort -n -o SRR7696207.namesort.bam SRR7696207.bam
$ samtools fixmate -m SRR7696207.namesort.bam SRR7696207.fixmate.bam
$ samtools sort -o SRR7696207.positionsort.fixmat.bam SRR7696207.fixmate.bam
$ samtools markdup -r SRR7696207.positionsort.fixmat.bam SRR7696207.rm.bam

看下以上的文件结构和大小

├── [ 136]  1
├── [ 136]  2
├── [   0]  config
├── [ 42K]  markdup.bam
├── [ 22G]  SRR7696207.bam
├── [5.4G]  SRR7696207.fixmate.bam
├── [5.2G]  SRR7696207.namesort.bam
├── [4.0G]  SRR7696207.positionsort.fixmat.bam
├── [3.8G]  SRR7696207.rm.bam
├── [2.6G]  SRR8517853.bam
├── [3.4G]  SRR8517854.bam
└── [6.9G]  SRR8517856.bam

最后查看下rm.bam

samtools view SRR7696207.rm.bam |less -S
SRR7696207.27107225     2193    chr1    10024   0       86M64H  chr5    10324   0       CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCT
SRR7696207.27728977     99      chr1    10025   17      91M18S  =       10025   91      TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTA
SRR7696207.27728977     147     chr1    10025   17      91M18S  =       10025   -91     TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTA
SRR7696207.4339191      163     chr1    10027   0       69M     =       10039   123     ACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC
SRR7696207.12487651     163     chr1    10028   0       86M     =       10028   81      CCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC
SRR7696207.9370306      83      chr1    10028   0       116M    =       10016   -128    CCCTAACCCTAACCCTAAACCTAACCCTAACCCTAACC
SRR7696207.12487651     83      chr1    10028   0       69S81M  =       10028   -81     TTTCTAGCAGTGGACTGCATACGTGTTCGCATACTGTG
SRR7696207.18904916     99      chr1    10030   0       81M69S  =       10030   79      CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCT
SRR7696207.18904916     147     chr1    10030   0       79M     =       10030   -79     CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCT
SRR7696207.4339191      83      chr1    10039   0       111M7S  =       10027   -123    ACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC
SRR7696207.20389847     99      chr1    10040   0       74M     =       10040   74      CCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC
SRR7696207.20389847     147     chr1    10040   0       76S74M  =       10040   -74     CATATTTTTGTTTTTTTTAATGTTACGGCGACCACCGA

看rm后的是否有差异

samtools flagstat SRR7696207.bam
samtools flagstat SRR7696207.rm.bam

结果如下

$ samtools flagstat SRR7696207.bam
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)
$ samtools flagstat SRR7696207.rm.bam
52114377 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
372636 + 0 supplementary
0 + 0 duplicates
52089646 + 0 mapped (99.95% : N/A)
51741741 + 0 paired in sequencing
25868532 + 0 read1
25873209 + 0 read2
51246880 + 0 properly paired (99.04% : N/A)
51699203 + 0 with itself and mate mapped
17807 + 0 singletons (0.03% : N/A)
319884 + 0 with mate mapped to a different chr
242709 + 0 with mate mapped to a different chr (mapQ>=5)

可以再试试-S参数

samtools markdup -S SRR7696207.sam SRR7696207.mk.bam

补充:samtoolsmarkdup操作的正确顺序

  1. The first sort can be omitted if the file is already name ordered
samtools sort -n -o namesort.bam example.bam
  1. Add ms and MC tags for markdup to use later
samtools fixmate -m namesort.bam fixmate.bam
  1. Markdup needs position order
samtools sort -o positionsort.bam fixmate.bam
  1. Finally mark duplicates
samtools markdup positionsort.bam markdup.bam

以上是简单的找变异流程,并不完善。下面是GATK找变异。

(0)

相关推荐

  • 软件介绍之Samtools

    咱们<生信技能树>的B站有一个lncRNA数据分析实战,缺乏配套笔记,所以我们安排了100个lncRNA组装案例文献分享,以及这个流程会用到的100个软件的实战笔记教程! 下面是100个l ...

  • 【直播】我的基因组25:用bcftools来call variation

    在这里,我们可以考虑比较对3种不同的bam文件来分别call variation的差异,探索对bam文件的不同过滤模式对snp calling的影响,分别是,原始比对的bam文件,去除低质量和多比对还 ...

  • 【直播】我的基因组19:根据比对结果来统计测序深度和覆盖度

    看来本次直播我的基因组分析流程效果还不错,不少朋友跟着自己动手开始分析全基因组测试数据了,值得表扬.其中有好几个朋友留言向我反映公司给的统计报告里面的覆盖度的问题,如下图: 我在第 8讲中写道,每条染 ...

  • sam转为bam文件报错

    怕什么报错呢,重来一次不就好了吗! 学员群有人提问,他们上完了我们的转录组课程,自己拿服务器去跑一个文献数据,有一个样本的sam转为bam文件报错,得到文件如下: 213M 2月  20 20:59  ...

  • 教程 | 简单粗暴的叶绿体基因组 SNP Calling 流程

    写在前面 最近主要忙一些植物群体基因组数据的项目.前面提过,赶时间,全基因组的 SNP Calling 使用 GATK 流程,还是需要跑上两三天.但这个还是耗时,并不一定能够赶上工期.于是我将目标转移 ...

  • 最新版针对RNA-seq数据的GATK找变异流程

    RNA-seq标准分析,我们已经讲解的太多了,表达矩阵到差异分析等下游生物学注释都没有啥新颖之处,融合基因和可变剪切算是出彩的地方,如果加上GATK找变异流程就更棒了,反正都使用了star软件进行序列 ...

  • 哪些人易被传染禽流感?中国科研团队在汉族人群中找出相关基因变异

    自然界中存在许多病毒,事实上它们大部分可以与人类相安无事.但有一些病毒不但能感染人类,还严重危害人类的健康.这其中,除了如今我们每天都在防范的新冠病毒外,一些流感病毒也不容忽视.为此,公共卫生专家和科 ...

  • 基因变异筛查预测乳腺癌风险

    A massive study of nearly 4,000 variants in a gene associated with cancer could help to pinpoint pe ...

  • 核辐射导致的基因变异,会不会遗传给下一代?

    本文来自微信公众号:学术经纬(ID:Global_Academia),作者:药明康德内容团队,原标题为<35年后再访切尔诺贝利,核辐射会如何影响基因变异?> 从广岛和长崎的两颗原子弹,到上 ...

  • 『艾滋病』基因变异:增加艾滋病毒感染 减慢艾滋病发病

    路透社纽约健康消息-根据<艾滋病>杂志电子网络版的一篇报导,人体基因的每个细微变异,都可引起感染艾滋病病毒的可能性增加:但一旦被感染,发展为艾滋病的过程却又可因此而被减慢. 美国国家过敏和 ...

  • 地球上基因变异人的“阿凡达人”

    人类的基因有时只要稍有突变,就会拥有一些神奇的能力.比如5号染色体上FAM134B基因的变异能让人们对疼痛免疫,变成不惧疼痛的超人:2号染色体上几个基因的突变可能会赋予你一双犹如戴着美瞳般眸色不同的眼 ...

  • 超越Pemigatinib?ARG087数据公布,针对FGFR2基因变异胆管癌

    Basilea Pharmaceutica公司于近日宣布,在在预先计划的II期试验的中期分析中,Derazantinib (ARG 087)证明了其抗肿瘤疗效,并达到了预先设定的阈值,可以按照计划对无 ...

  • 一种可抵抗疟疾的基因变异

    <正> 澳大利亚墨尔本Walter andEliza Hall医学研究所(WEHI)的AlanF.Cowman博士及其同事在2003年1月的<自然医学>杂志上发表文章称,他们在 ...

  • 美国发现一基因变异猩猩,手指已进化成人类一样,它能否持续进化

    2019年,美国亚特兰大动物园一工作人员在自己的社交平台上分享了一些照片,里面是一只名叫Anaka的黑猩猩,这天是它6岁生日,工作人员们正在给它过生日,非常喜庆,然而网络却因为这些照片"炸锅 ...