最近在學(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 指定輸出文件的位置及文件名
這個(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í)間再來看這塊吧)