Samtools的一个小实例

最近在学习利用转录组数据进行变异检测的相关分析方法,找到了一篇关于samtools的教程SAMtools - Primer / Tutorial,原文提供分析过程用到的数据,非常好的学习素材,简单记录自己的重复过程。(前两天试着用conda安装了samtools-1.8-3, 但是使用的时候遇到了报错

暂时先不管它了,用之前安装的samtools-1.5来重复本次教程;conda之前好长一段时间都不能用,这几天又突然可以用了,搞不懂什么原因)

PS:看完原文教程以后发现内容非常熟悉,仔细想想原来是老师上课讲过的一个实例!

典型的二代测序数据分析流程主要包括五个步骤(A typical work-flow for next-generation sequence analysis is usually composed of the following steps:)

1 DNA is extracted from a sample(DNA提取)

2 DNA is sequenced.(测序)

3 Raw sequencing reads are aligned to a reference genome(reads比对到参考基因组)

4 Aligned reads are evaluated and visualized.(比对结果评估和可视化)

5 Genomic variants, variants, including single nucleotide polymorphisms (SNPs), small insertions and deletions are identified(SNP,Indel检测)

步骤4,5用到 samtools

本次教程的流程(workflow)主要包括5个步骤

1 Align the reads to the reference E.coli genome with Bowtie2(比对reads到参考基因组)

2 convert the aligned reads from the SAM file format to BAM(将SAM格式转换为BAM格式)

3 sort and index the BAM file (index该怎么翻译呢?)

4 identify genomic variants(变异检测)

5 visualize the reads and genomic variants(可视化)

(接下来原文用wgsim模拟生成了原始的reads,直接跳过这一部分)

Bowtie2将reads比对到参考基因组之前,首先要index参考基因组

Bowtie 2: Manual 这个链接是Bowtie2的帮助文档,内容非常详细,只是要看懂还得费些力气

index参考基因组用到的命令是 bowtie2-bulid (命令后的参数很多自己也还有很多参数没有搞懂),通常直接用默认参数,用到的命令为 bowtie2-build input_file.fasta output_index_name


试了好几次一直报错也没发现问题出在哪,最后才知道实验室的服务器没有存储空间了,甚至用mkdir命令新建目录都告诉我permission denied了

接下来是用bowtie2命令将reads比对到参考基因组

-x参数指定index后的参考基因组位置 -U reads所在的位置 -S 指定输出文件的位置及文件名

比对结果,应该是aligned exactly 1 time的值越大越好

这个是单端数据,如果是PE 用-1和-2分别指定reads

接下来原文介绍了sam文件的格式,记得生物信息学课程考试考了这个题,自己好像没有写出来(这个格式现在的自己也没有太明白,得抽时间仔细看好好理解啦!)

我们的目标是鉴定基因组变异,第一步要做的是将sam转换为bam,下游的分析需要BAM文件作为输入(our goal is to identity the set of genomic variants within the E.coli data set. To do so, our first step is to convert the SAM file to BAM. This is an important prerequisite, as all the downstream steps, including the identification of genomic variants visualization of reads, require BAM as input)

转换sam为bam用到samtools view 命令

接下来是sort和index bam文件,分别用到samtools sort 和 samtools index

sort的时候可能因为samtools的版本不一样重复原文命令会遇到报错,需要我们自己加上一个 -o 参数指定输出文件

现在有了bam文件和参考基因组,试着将比对结果可视化一下,用到IGV


通过2导入参考基因组fasta文件,通过1导入index后的bam文件(index时生成的后缀为fai的文件也必须同时存在),通过4和5 可以缩小和放大参考基因组显示的区域,通过改变3破折号前后的数字可以选择要展示的基因组区域(IGV好像还有很多小工具有机会可以探索一下)

IGV查看比对结果所用到的文件  百度云链接密码 il9q

第一次接触IGV是在跟着生信技能树论坛学习转录组数据分析的时候,当时用IGV查看比对的结果没有成功,后来也没有在尝试,回头想一想,都过去一年多的时间了......

变异检测

第一步使用samtools mpileup 计算 the genotype likelihoods supported by the aligned reads in our sample

然后用bcftools来检测变异(另外抽时间再来看这块吧)

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,530评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 86,403评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,120评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,770评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,758评论 5 367
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,649评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,021评论 3 398
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,675评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,931评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,659评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,751评论 1 330
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,410评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,004评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,969评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,042评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,493评论 2 343

推荐阅读更多精彩内容