宏基因組實戰7. bwa序列比對, samtools查看, bedtools豐度統計

2021-01-14 宏基因組
前情提要

如果您在學習本教程中存在困難，可能因為缺少背景知識，建議先閱讀本系統前期文章

測試數據

劉博士幫助把測試數據建立了一個百度雲同步共享文件夾，有非常多的好處，請讀完下文再決定是否下載：

下載被牆的數據；很多數據存在google, amazon的部分伺服器國內無法直接下載，而伺服器一般科學上網不方便，下載數據困難。大家下載失敗的數據請到共享目錄中查找；

預下載好的軟體、資料庫；有很多需要下載安裝、註冊的軟體(在線安裝包除外)，其實已經在共享目錄了，節約小夥伴申請、下載的時間；

數據同步更新；任何筆記或教程不可避免的有些錯誤、或不完善的地方，後期通過大家的測試反饋問題，我可以對教程進行改進。共享目錄不建議全部下載或轉存，因為文件體積非常大(目前>18G)，而且還會更新。你轉存的只是當前版本的一個備份，就不會再更新了。建議直接在連結中每次逐個下載需要的文件，也對文件有一個認識過程。

方便結果預覽和跳過問題步驟；伺服器Linux在不同平臺和版本下，軟體安裝和兼容性問題還是很多的，而且用戶的權限和經驗也會導致某些步驟相關軟體無法成功安裝(有問題建議選google、再找管理員幫助；想在群裡提問或聯繫作者務必閱讀《如何優雅的提問》)。在百度雲共享目錄中，有每一步的運行結果，讀者可以下載查看分析結果，並可基於此結果進一步分析。不要糾結於某一步無法通過，重點是了解整個流程的分析思路。

最後送上本教程使用到的所有文件同步共享文件夾連結：http://pan.baidu.com/s/1hsIjosk 密碼：y0tb 。

bwa序列比對

安裝bwa

# 工作目錄，根據個人情況修改wd=~/test/metagenome17cd $wdsudo apt-get install bwa samtools

下載數據

mkdir mappingcd mapping# 可以接之前的學習，或在百度雲中下載cp ../data/*.pe.fq.gz ./ # 引處如果ln硬鏈解壓不允許# 逐個解壓，直接gunzip *.gz批量也不成功for file in *fq.gzdogunzip $filedone# 無法下載找百度雲curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/subset_assembly.fa.gzgunzip subset_assembly.fa

序列比對

# 建索引bwa index subset_assembly.fa# 雙端合併序列，使用-p參數比對for i in *fqdo# unrecognized command 'mem' 版本過舊，qiime索引了舊版本，bwa改為絕對路/usr/bin/bwa mem -p subset_assembly.fa $i > ${i}.aln.sam done

samtools操作比對結果

samtools可視化比對結果

# 參考序列建索引samtools faidx subset_assembly.fa# 壓縮sam為bam用於可視化for i in *.samdo samtools import subset_assembly.fa $i $i.bam samtools sort $i.bam -o $i.bam.sorted.bam samtools index $i.bam.sorted.bamdone# 可視化# 按contig的reads數量排序，找高豐度的查看grep -v ^@ SRR1976948.abundtrim.subset.pe.fq.aln.sam | cut -f 3 | sort | uniq -c | sort -n | tail# Pick one e.g. k99_13588.# 查看k99_13588序列400bp開始samtools tview SRR1976948.abundtrim.subset.pe.fq.aln.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:400# 方向可以上下左右移動查看，q退出# 另一個窗口打開另一個樣品，便於比較samtools tview SRR1977249.abundtrim.subset.pe.fq.aln.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:400

在兩個終端(Xshell)中用samtools查看不同樣品的比對結果序列分析情況

bedtools豐度估計

我們將會用比對結果估計組裝序列的覆蓋度

使用bedtools比對，常用 http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html

# 安裝程序sudo apt-get install bedtools# 使用genomeCoverageBed定量基因for i in *sorted.bam    do    genomeCoverageBed -ibam $i > ${i/.pe*/}.histogram.tabdone

我們看一下結果：

k99_5    6    87    258    0.337209k99_5    7    18    258    0.0697674k99_5    8    11    258    0.0426357

Contig name

Depth of coverage 覆蓋深度

Number of bases on contig depth equal to column 2

Size of contig (or entire genome) in base pairs

Fraction of bases on contig (or entire genome) with depth equal to column 2

To get an esimate of mean coverage for a contig we sum (Depth of coverage) * (Number of bases on contig) / (Length of the contig). We have a quick script that will do this calculation.

計算平均深度

wget https://raw.githubusercontent.com/ngs-docs/2017-cicese-metagenomics/master/files/calculate-contig-coverage.py# 安裝過可跳過sudo pip install pandasfor hist in *histogram.tabdo    python calculate-contig-coverage.py $histdone

產生結果文件 SRR1976948.abundtrim.subset.histogram.tab.coverage.tab：

k99_100    18.200147167k99_1000    52.6492890995k99_10000    4.97881916865k99_10008    16.7718223583k99_1001    62.2604006163

可選分析：末質控數據結果如何？

# 下載原始數據curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gzcurl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz# 如果有，可複製過來cp ../data/SRR1976948_?.fastq.gz .# 解壓變只提取前20萬行gunzip -c SRR1976948_1.fastq.gz | head -800000 > SRR1976948.1gunzip -c SRR1976948_2.fastq.gz | head -800000 > SRR1976948.2# 比對bwa aln subset_assembly.fa SRR1976948.1 > SRR1976948_1.untrimmed.saibwa aln subset_assembly.fa SRR1976948.2 > SRR1976948_2.untrimmed.saibwa sampe subset_assembly.fa SRR1976948_1.untrimmed.sai SRR1976948_2.untrimmed.sai SRR1976948.1 SRR1976948.2 > SRR1976948.untrimmed.sam# 壓縮、排序、索引i=SRR1976948.untrimmed.samsamtools import subset_assembly.fa $i $i.bamsamtools sort $i.bam -o $i.bam.sorted.bamsamtools index $i.bam.sorted.bam# 查看samtools tview SRR1976948.untrimmed.sam.bam.sorted.bam subset_assembly.fa -p k99_13588:500

比較這些沒有trim質控的數據，看看和原來的結果有什麼不同？

猜你喜歡寫在後面

為鼓勵讀者交流、快速解決科研困難，我們建立了「宏基因組」專業討論群，目前己有國內外七十多位PI，七百多名一線科研人員加入。參與討論，獲得專業指導、問題解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註「姓名-單位-研究方向-職務」。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍末解決群內討論，問題不私聊，幫助同行。

學習16S擴增子、宏基因組科研思路和分析實戰，關注「宏基因組」

相關焦點

  • Bedtools使用簡介
    提取給定位置的FASTA序列    maskfasta     Use intervals to mask sequences from a FASTA file.    @ehbio:~/bedtools$ curl -O ${url}/gwas.bedct@ehbio:~/bedtools$ curl -O ${url}/genome.txtct@ehbio:~/bedtools$ curl -O ${url}/hesc.chromHmm.bed交集 (intersect)查看輸入文件，bed格式，至少三列，分別是染色體，起始位置(
  • 如何進行基因組序列比對?
    VN:0.7.13-r1126   CL: bwa mem hg19.fa read1.fq  read2.fqID: 軟體名稱VN: 軟體版本號CL: 命令行 @HD用SO來記錄比對後的reads其排列順序示例：SO:unsortedSO:queryname
  • 處理SAM、BAM你需要Samtools
    sam是短序列比對默認的標準格式，由sanger制定，是以TAB為分割符的文本格式。主要應用於測序序列mapping到基因組上的結果表示，另外也可以表示其他的多重比對結果。，但在基因組中一般作為結構變異的證據【下面👇的圖1是原始數據(其中Coor是坐標的簡寫；ref是參考序列)；圖2是使用bwa、bowtie2、hisat2等軟體比對的結果sam文件】【其中r001/1與r001/2是paired end read； r003是嵌合read，有兩個記錄】
  • nanopore測序技術專題(十九):利用samtools處理sam格式文件
    第四列：為在參考序列上的位置第五列：比對的質量值，MAPQ第六列：代表比對結果的CIGAR字符串第七列：mate比對到的染色體號，若是沒有mate，則是*第八列：比對到參考序列上的第一個鹼基位置第九列：Template的長度，第十列：為read的序列第十一列：為ASCII碼格式的序列質量；利用samtools處理sam/bam
  • samtools命令大全
    samtools是一個用於操作sam和bam文件（通常是短序列比對工具如bwa,bowtie2,hisat2,tophat2等等產生的，
  • GATK4.0全基因組和全外顯子組分析實戰(修改版)
    index -a bwtsw -p hg19 hg19.fa  &程序運行時間較長，建議使用nohup 或者 screen最終生成文件  hg19.amb  hg19.ann  hg19.bwt  hg19.pac hg19.sabwa比對bwa mem -t 4 -M -R 『
  • 一個全基因組重測序分析實戰
    /bio-bwa/files/bwa-0.7.15.tar.bz2tar xvfj bwa-0.7.15.tar.bz2 # x extracts, v is verbose (details of what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 files
  • 生信人必會之samtools
    很多情況下需要有bai文件的存在，特別是顯示序列比對情況下。比如samtool的tview命令就需要；gbrowse2顯示reads的比對圖形的時候也需要。鍵來查看。PS:建議不要學這個命令了，還是直接用IGV看吧！扔給你BAM比對結果！！
  • 多快好省的宏基因組研究技巧
    下圖是宏基因組測序數據中比對到人類基因組的序列比例，根據樣本類型不同而不同。HUMAnN2，其首先使用MetaPhlAn2進行物種分類（關於這個軟體我們在前面物種組成部分已經講過），並提取相應物種參考基因組用於比對，未比對上的用於進一步和uniref資料庫進行蛋白質序列比對。
  • 快速看懂腸道菌群宏基因組測序分析報告
    -K55)，其他組裝參數與單樣品組裝參數相同；6) 將混合組裝的Scaffolds從N連接處打斷，得到不含N的Scaftigs序列；7) 對於單樣品和混合組裝生成的Scaftigs，過濾掉500bp以下的片段，並進行統計分析和後續基因預測；各樣品組裝結果基本信息統計說明：SampleID為樣品名稱；Total Length表示組裝得 Scaftigs
  • Bioinformatic Data Skills 學習專題(7) 比對數據的小實操
    找特定區域上的比對情況作者再github項目上分享了1000 Genomes Project Consortium的數據來實戰，具體可以到https://github.com/vsbuffalo/bds-files/tree/master/chapter-11-alignment 下載samtools index NA12891_CEU_sample.bam
  • 宏基因組序列物種分類之kraken 1/2和Bracken的使用
    人的轉錄組測序完，想快速看看人、微生物的序列的比例?元/宏基因組測序完，想快速獲得樣本中物種的豐度信息？Kraken 1Kraken 1是2014年Wood DE在Genome Biology發表的宏基因組序列分類軟體，能夠快速對宏基因樣品中的reads進行分類。Kraken在序列比對環節基於精確k-mer匹配和精簡資料庫的方法，採取精確匹配，其核心是Kraken有一種特殊資料庫，用以預先計算序列中包含的特殊的Kmer序列。
  • 常用在線序列比對工具
    序列比對算法。但是他大爺畢竟是他大爺，如果想比對兩條序列全局情況，Needle仍是不錯的選擇，至於Smith-Waterman算法，在許多二代比對算法裡仍可見，比如bwa-sw算法，Minimap2計算overlap時使用的也是Smith-Waterman算法；Smith-Waterman也擁有CUDA版本程序。
  • 生信第二步:建立索引及比對
    所以mRNA反轉的cDNA如果比對不到參考序列，會被分開，重新比對一次，判斷中間是否有內含子。連結：https://www.jianshu.com/p/681e02e7f9af第三個問題：如何比對,有哪些軟體的選擇？
  • Nature子刊:HUMAnN2實現宏基因組和宏轉錄組種水平功能組成分析
    Gregory Caporaso,  Nicola Segata簡介HUMAnN是一款快速宏基因組功能組成定量工具，第一版2012年發表於PLoS computational biology (當年還沒有影響因子，最新17年3.9，歷史最高14年4.6)，截止2018年11月1號Google Scholar統計引用557次。