Samtools的一個(gè)小實(shí)例

最近在學(xué)習(xí)利用轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行變異檢測的相關(guān)分析方法尺碰,找到了一篇關(guān)于samtools的教程SAMtools - Primer / Tutorial蹈矮,原文提供分析過程用到的數(shù)據(jù)注服,非常好的學(xué)習(xí)素材舆驶,簡單記錄自己的重復(fù)過程拓轻。(前兩天試著用conda安裝了samtools-1.8-3, 但是使用的時(shí)候遇到了報(bào)錯(cuò)

暫時(shí)先不管它了,用之前安裝的samtools-1.5來重復(fù)本次教程批幌;conda之前好長一段時(shí)間都不能用础锐,這幾天又突然可以用了,搞不懂什么原因)

PS:看完原文教程以后發(fā)現(xiàn)內(nèi)容非常熟悉荧缘,仔細(xì)想想原來是老師上課講過的一個(gè)實(shí)例皆警!

典型的二代測序數(shù)據(jù)分析流程主要包括五個(gè)步驟(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.(比對結(jié)果評估和可視化)

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

步驟4,5用到 samtools

本次教程的流程(workflow)主要包括5個(gè)步驟

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格式轉(zhuǎn)換為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 這個(gè)鏈接是Bowtie2的幫助文檔绸罗,內(nèi)容非常詳細(xì)意推,只是要看懂還得費(fèi)些力氣

index參考基因組用到的命令是 bowtie2-bulid (命令后的參數(shù)很多自己也還有很多參數(shù)沒有搞懂),通常直接用默認(rèn)參數(shù)珊蟀,用到的命令為 bowtie2-build input_file.fasta output_index_name


試了好幾次一直報(bào)錯(cuò)也沒發(fā)現(xiàn)問題出在哪菊值,最后才知道實(shí)驗(yàn)室的服務(wù)器沒有存儲空間了,甚至用mkdir命令新建目錄都告訴我permission denied了

接下來是用bowtie2命令將reads比對到參考基因組

-x參數(shù)指定index后的參考基因組位置 -U reads所在的位置 -S 指定輸出文件的位置及文件名

比對結(jié)果,應(yīng)該是aligned exactly 1 time的值越大越好

這個(gè)是單端數(shù)據(jù)俊性,如果是PE 用-1和-2分別指定reads

接下來原文介紹了sam文件的格式略步,記得生物信息學(xué)課程考試考了這個(gè)題描扯,自己好像沒有寫出來(這個(gè)格式現(xiàn)在的自己也沒有太明白定页,得抽時(shí)間仔細(xì)看好好理解啦!)

我們的目標(biāo)是鑒定基因組變異绽诚,第一步要做的是將sam轉(zhuǎn)換為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)

轉(zhuǎn)換sam為bam用到samtools view 命令

接下來是sort和index bam文件典徊,分別用到samtools sort 和 samtools index

sort的時(shí)候可能因?yàn)閟amtools的版本不一樣重復(fù)原文命令會(huì)遇到報(bào)錯(cuò),需要我們自己加上一個(gè) -o 參數(shù)指定輸出文件

現(xiàn)在有了bam文件和參考基因組恩够,試著將比對結(jié)果可視化一下卒落,用到IGV


通過2導(dǎo)入?yún)⒖蓟蚪Mfasta文件,通過1導(dǎo)入index后的bam文件(index時(shí)生成的后綴為fai的文件也必須同時(shí)存在)蜂桶,通過4和5 可以縮小和放大參考基因組顯示的區(qū)域儡毕,通過改變3破折號前后的數(shù)字可以選擇要展示的基因組區(qū)域(IGV好像還有很多小工具有機(jī)會(huì)可以探索一下)

IGV查看比對結(jié)果所用到的文件? 百度云鏈接密碼 il9q

第一次接觸IGV是在跟著生信技能樹論壇學(xué)習(xí)轉(zhuǎn)錄組數(shù)據(jù)分析的時(shí)候,當(dāng)時(shí)用IGV查看比對的結(jié)果沒有成功扑媚,后來也沒有在嘗試腰湾,回頭想一想,都過去一年多的時(shí)間了......

變異檢測

第一步使用samtools mpileup 計(jì)算 the genotype likelihoods supported by the aligned reads in our sample

然后用bcftools來檢測變異(另外抽時(shí)間再來看這塊吧)

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末疆股,一起剝皮案震驚了整個(gè)濱河市费坊,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌旬痹,老刑警劉巖附井,帶你破解...
    沈念sama閱讀 222,681評論 6 517
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異两残,居然都是意外死亡永毅,警方通過查閱死者的電腦和手機(jī),發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,205評論 3 399
  • 文/潘曉璐 我一進(jìn)店門人弓,熙熙樓的掌柜王于貴愁眉苦臉地迎上來沼死,“玉大人,你說我怎么就攤上這事票从÷瘢” “怎么了?”我有些...
    開封第一講書人閱讀 169,421評論 0 362
  • 文/不壞的土叔 我叫張陵峰鄙,是天一觀的道長浸间。 經(jīng)常有香客問我,道長吟榴,這世上最難降的妖魔是什么魁蒜? 我笑而不...
    開封第一講書人閱讀 60,114評論 1 300
  • 正文 為了忘掉前任,我火速辦了婚禮,結(jié)果婚禮上兜看,老公的妹妹穿的比我還像新娘锥咸。我一直安慰自己,他們只是感情好细移,可當(dāng)我...
    茶點(diǎn)故事閱讀 69,116評論 6 398
  • 文/花漫 我一把揭開白布搏予。 她就那樣靜靜地躺著,像睡著了一般弧轧。 火紅的嫁衣襯著肌膚如雪雪侥。 梳的紋絲不亂的頭發(fā)上,一...
    開封第一講書人閱讀 52,713評論 1 312
  • 那天精绎,我揣著相機(jī)與錄音速缨,去河邊找鬼。 笑死代乃,一個(gè)胖子當(dāng)著我的面吹牛旬牲,可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播搁吓,決...
    沈念sama閱讀 41,170評論 3 422
  • 文/蒼蘭香墨 我猛地睜開眼原茅,長吁一口氣:“原來是場噩夢啊……” “哼!你這毒婦竟也來了擎浴?” 一聲冷哼從身側(cè)響起员咽,我...
    開封第一講書人閱讀 40,116評論 0 277
  • 序言:老撾萬榮一對情侶失蹤,失蹤者是張志新(化名)和其女友劉穎贮预,沒想到半個(gè)月后贝室,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 46,651評論 1 320
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡仿吞,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 38,714評論 3 342
  • 正文 我和宋清朗相戀三年滑频,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片唤冈。...
    茶點(diǎn)故事閱讀 40,865評論 1 353
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡峡迷,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出你虹,到底是詐尸還是另有隱情绘搞,我是刑警寧澤,帶...
    沈念sama閱讀 36,527評論 5 351
  • 正文 年R本政府宣布傅物,位于F島的核電站夯辖,受9級特大地震影響,放射性物質(zhì)發(fā)生泄漏董饰。R本人自食惡果不足惜蒿褂,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 42,211評論 3 336
  • 文/蒙蒙 一圆米、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧啄栓,春花似錦娄帖、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,699評論 0 25
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至桂肌,卻和暖如春数焊,著一層夾襖步出監(jiān)牢的瞬間,已是汗流浹背崎场。 一陣腳步聲響...
    開封第一講書人閱讀 33,814評論 1 274
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留遂蛀,地道東北人谭跨。 一個(gè)月前我還...
    沈念sama閱讀 49,299評論 3 379
  • 正文 我出身青樓,卻偏偏與公主長得像李滴,于是被迫代替她去往敵國和親螃宙。 傳聞我的和親對象是個(gè)殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,870評論 2 361

推薦閱讀更多精彩內(nèi)容

  • 五點(diǎn)打坐1小時(shí) 6點(diǎn)哈他1小時(shí)所坯, 7點(diǎn)早餐:??雞蛋 8點(diǎn)讀經(jīng)書1小時(shí) 9點(diǎn)瑜伽視頻1小時(shí) 10點(diǎn)看1小時(shí)瑜伽書 1...
    A快樂我做主閱讀 202評論 0 0
  • 請點(diǎn)擊此處輸入圖片描述 Q:當(dāng)今社會(huì)芹助,為什么越來越多的女人主動(dòng)提出離婚呢堂湖? A:女人主動(dòng)提出離婚,做老公的根本就沒...
    飄雨桐V閱讀 793評論 0 0
  • ScrollView 一個(gè)包裝了平臺的ScrollView(滾動(dòng)視圖)的組件状土,同時(shí)還集成了觸摸鎖定的“響應(yīng)者”系統(tǒng)...
    一畝水塘閱讀 2,806評論 3 11