重測序之BSA

劉小澤寫于18.7.18
表型與基因型相關(guān)聯(lián)（基因定位）是遺傳學(xué)上的重要問題。

什么是重測序Resequencing？

我們使用的基因組是怎么來的？通過Denovo測序，也就是從頭測序，拼接得到基因組

那么重測序呢？已有的基因組發(fā)表出來以后，重新進行測序就不需要拼接，直接可以對測序reads進行短序列比對。簡而言之就是對基因組序列已知的物種個體進行基因組測序。本質(zhì)就是找變異。

通過這種方法，可以比較個體或者群體差異，建立該物種的變異數(shù)據(jù)庫，挖掘重要性狀的關(guān)鍵候選基因。利用重測序可以進行變異檢測(SNP Calling)、遺傳圖譜構(gòu)建、群體進化以及全基因組關(guān)聯(lián)分析（GWAS）。

由于基因組的測序覆蓋度不同，可分為全基因組重測序和簡化基因組重測序。全基因組重測序（WGS）是在全基因組水平上檢測動植物重要性狀相關(guān)的變異位點，可以檢測SNP(單核苷酸多態(tài)性)、InDel（小片段插入缺失）、SV（結(jié)構(gòu)變異）、CNV（拷貝數(shù)變異）、轉(zhuǎn)座子變異等。

利用鳥槍法建立350bp文庫，測序深度10-30X，適用于有參物種，成本較高(SNP/InDel >10X; SV > 20X; CNV > 30X, 轉(zhuǎn)座子> 20X)增加測序深度可以提高基因組覆蓋度和變異檢出率

簡化基因組測序（GBS, Genotyping by Sequencing） 將基因組DNA進行酶切，然后對酶切片段序列進行測序，一般能捕獲全基因組的1%序列，主要用于檢測SNP，成本低適合大規(guī)模樣品標記篩選。

什么是BSA？

BSA是一種檢測位于特定染色體片段上標記的方法，又叫Mapping-by-sequencing。理論上，任何一對具有相對形狀的親本雜交產(chǎn)生的分離后代都適用，也就是針對質(zhì)量性狀單基因或數(shù)量性狀主效基因的定位。但BSA一次只能分析一個性狀

建庫時兩個親本分別建庫，子代分別混池建庫，共四個，文庫片段大小350bp測序深度一般親本>10X，子代混池>20X

一般要經(jīng)過以下幾步：

  1. 發(fā)現(xiàn)相對性狀：可能是人為誘變，可能是自然發(fā)現(xiàn)

  2. 構(gòu)建群體: 異交(outcross)或回交(backcross)

  3. 選取目標性狀群體

  4. 測序分析

如何分析？

拿回來數(shù)據(jù)，第一步總是要對數(shù)據(jù)過濾清洗，去接頭去不合格reads；第二步是序列比對，常用bwa、bowtie2；第三步SNP calling， GATK(最全最好用的)、samtools+bcftools（易上手的）；第四步作圖確定可能區(qū)域；第五步變異位點注釋

使用的軟件：

主要是SHOREmap 、samtools、bcftools、fastqc、multiqc、bwa，可以先把軟件安裝好，最快的方法是conda install -p /your conda path/ tools 。

下載測試數(shù)據(jù)：

這一次分析擬南芥回交BC群體數(shù)據(jù)以及參考基因組、注釋文件

mkdir -p ~/bwa/data && cd ~/bwa/data
axel http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads1.fq.gz
axel http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads2.fq.gz
axel http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads1.fq.gz
axel http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads2.fq.gz
# reference genome & gff in dir reference 下載到reference下
mkdir -p ~/bwa/reference && cd ~/bwa/reference
wget -c https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas
wget -c https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff
?
#chr_sizes
wget http://bioinfo.mpipz.mpg.de/shoremap/data/software/chrSizes.txt</pre>

分析過程：

  1. 質(zhì)控

    進入到數(shù)據(jù)文件夾，例如data下

    ls *.fq | xargs -i fastqc {}
multiqc .
1.png
2.png

這應(yīng)該是清洗以后的數(shù)據(jù)了，質(zhì)量還不錯。如果自己拿回來的數(shù)據(jù)，有時需要去接頭，去掉前幾bp的低質(zhì)量序列

  1. 序列比對

    重測序中的比對使用最多就是bwa，分為兩步：先對參考基因組構(gòu)建索引index，然后使用bwa mem比對；【當(dāng)然構(gòu)建索引使用bowtie2軟件的bowtie-build命令也是可以的】比對的結(jié)果是sam格式，我們這里用samtools sort一步將管道輸入給它的sam文件排序并轉(zhuǎn)換格式為bam文件【你可能注意到還有一些start、end、echo之類的命令，他們是為了提示用戶進行到了哪一步，并且計算運行時間；這些都是可有可無的，屬于點綴】

3.png

  1. SNP calling這是基于貝葉斯定理，你肯定眼熟～后期講統(tǒng)計概率先講它！太重要了！

    4.png

    假如我們之前比對的短序列和參考基因組有差異了，那么還不能確定他就是一個變異，作為候選人，需要接受多方考證：一個是測序準確度，測不準也有可能會有變異；一個是比對軟件的精度，都是人開發(fā)的，難免有誤差：一個是抽樣的誤差，就是群體選材時產(chǎn)生的；如何確定它是不是一個真實的變異呢？這里貝葉斯定理就派上了大用場～

    比如基因組上一個位點有50條reads比對上，但是其中的40條都與參考基因組不同，那么很有可能他是一個真的變異位點，至于什么變異下面再分析

    這里用的samtools是最新版本1.8，mplieup主要功能主要是生成BCF、VCF文件。這里-v指定輸出為vcf文件，否則默認生成pileup格式文件；-f指定構(gòu)建過索引的參考序列；-u是uncompressed，輸出未壓縮的VCF文件；

5.png

  1. SHOREmap確定候選區(qū)域

    SHOREmap的安裝還是比較復(fù)雜的，需要有root權(quán)限去安裝兩個動態(tài)庫首先，確定自己的系統(tǒng)是不是64位uname -m 顯示X86_64就是64位；其次，安裝依賴的庫DISLIN ，而doslin自己又需要兩個動態(tài)庫

#檢查安裝doslin的兩個動態(tài)庫
    #檢查libXm.so*是否安裝，在命令行輸入
    ls /usr/lib/libXm.so*
    #如果顯示'ls: cannot access /usr/lib/libXm.so*: No such file or directory‘，就需要安裝【接下來需要root權(quán)限】
    sudo apt-get update
    sudo apt-get install libmotif4
    #再檢查一遍
    ls /usr/lib/libXm.so*
    ?
    #檢查libXt.so是否安裝
    ls /usr/lib/x86_64-linux-gnu/libXt.so 
    #32位系統(tǒng)輸入：ls /usr/lib/i386-linux-gnu/libXt.so
    #出現(xiàn)'ls: cannot access /usr/lib/x86_64-linux-gnu/libXt.so: No such file or directory‘，就需要安裝
    sudo apt-get install libxt-dev
    #再檢查一遍
    #安裝完doslin需要的動態(tài)庫后，進入軟件目錄，開始下載doslin
    cd /path/to/software
    # 下載
    wget ftp://ftp.gwdg.de/pub/grafik/dislin/linux/i586_64/dislin-11.1.linux.i586_64.tar.gze
    # 解壓縮
    tar -zxvf dislin-11.0.linux.i586_64.tar.gz
    cd dislin-11.0
    # 加入系統(tǒng)路徑
    mkdir -p $HOME/biosoft/dislin
    DISLIN=$HOME/biosoft/dislin
    export DISLIN
    # 開始安裝
    ./INSTALL
    # 復(fù)制dislin_d.h 到dislin文件夾中
    cp ./examples/dislin_d.h $DISLIN


    然后才開始正題——安裝SHOREmap，當(dāng)前最新版本是3.6版本
    #里面的安裝路徑需要個性化自己定制哦
    cd /your software path/
    wget http://bioinfo.mpipz.mpg.de/shoremap/SHOREmap_v3.6.tar.gz
    # 替換SHOREmap下的dislin的一些文件
    tar -zxvf SHOREmap_v3,6
    rm dislin/*dislin_d.*
    cp $DISLIN/*dislin_d.* dislin
    # 編輯/etc/profile[root用戶]或.bashrc[非root]
    vi ~/.bashrc
    export LD_LIBRARY_PATH=/path/to/SHOREmap_v3.6/dislin
    # 退出保存.bashrc
    source ～/.bashrc
    # 到之前安裝的文件夾下
    cd & cd src/SHOREmap_v3.6
    #進行編譯(可選，如果沒有g(shù)++)
    #sudo apt-get install build-essential
    make
    #【編譯過程如果出錯，有可能動態(tài)庫文件路徑設(shè)置不正確，把它放到/etc/profile或者/.bashrc的上面】
    # 將編譯文件拷貝到習(xí)慣的文件夾中，然后添加執(zhí)行路徑
    cp SHOREmap ../../biosoft/SHOREmap_v3.6
    echo "export $HOME/bisoft/SHOREmap_v3.6" >> ~/.bashrc
    source ～/.bashrc
    # 檢查是否安裝成功
    SHOREmap
6.png

以下是官網(wǎng)給出的針對Outcross和Backcross群體的不同分析流程，但是都要求先轉(zhuǎn)換格式【我們這里的測試流程是針對的Backcross】

6_1.png
  • 轉(zhuǎn)換格式：將vcf轉(zhuǎn)換成軟件能認識的格式基本兩個參數(shù)就是--marker 和--folder【注意：因為這里是用的兩臺服務(wù)器，前面跑的那臺沒有root權(quán)限，又注冊了一個騰訊云服務(wù)器，獲得了權(quán)限后安裝的各種庫，接下來的分析也用的是云服務(wù)器，因此代碼可能和上面有出入】

SHOREmap convert --marker ./variant/BC_bg_raw_variant.vcf --folder background/
SHOREmap convert --marker ./variant/BC_fg_raw_variant.vcf --folder foreground/
#生成三個文件??</pre>

7.png

  • 提取候選分子標記的一致性信息:SHOREmap extract

    四個參數(shù)--chrsizes 、--folder、--marker 、--consen 分別對應(yīng)染色體長度信息、指定輸出目錄、上面轉(zhuǎn)換后的

    #background
    SHOREmap extract --chrsizes chrSizes.txt --folder ../background/ --marker ../background/1_converted_variant.txt --consen ../background/1_converted_consen.txt -verbose
    #foreground
    SHOREmap extract --chrsizes chrSizes.txt --folder ../foreground/ --marker ../foreground/1_converted_variant.txt --consen ../foreground/1_converted_consen.txt -verbose
    #生成 extracted_consensus_0.txt</pre>
x_8.png

  • 回交分析SHOREmap backcross用來分析回交作圖群體得到的重組后代混池數(shù)據(jù)，想過濾出所有參考基因組和測序池之間的不同部分用于找到突變點特異部分。如何確保不是自然變異或者是測序錯誤而產(chǎn)生的特異部分呢？要求測序池中選擇的部分要多次出現(xiàn)在親本或背景中。再根據(jù)前景/背景的堿基識別,base calls,質(zhì)量/覆蓋率/等位基因等信息，確定是否把保留的SNP位點作為分子標記

    必選的參數(shù)有三個：

    1. 染色體大小文件--chrsizes，包括兩行，一是染色體位置，二是大小；

    2. --marker通過SHOREmap convert生成的候選marker文件 converted_variant.txt

      x_9.png

      這幾列分別代表：項目名稱、染色體ID、SNP標記位點位置、參考序列上的堿基、變異的堿基、可變堿基質(zhì)量（0-40）【也是用Phred換算】、符合這種堿基改變的reads數(shù)、前一項占總覆蓋度的百分比、覆蓋到SNP位置的比對平均命中數(shù)（也就是這個區(qū)間能檢測到SNP的位點數(shù)）

    3. --folder 輸出文件夾

SHOREmap backcross --chrsizes chrSizes.txt --marker ../foreground/1_converted_variant.txt --consen extracted_consensus_0.txt --folder ./ -plot-bc --marker-score 40 --marker-freq 0.0 --min-coverage 10 --max-coverage 80 -bg ../background/1_converted_variant.txt --bg-cov 1 --bg-freq 0.0 --bg-score 1 --non-EMS --cluster 1 --marker-hit 1 -verbose

顯示的結(jié)果及文件：
x_10.png

結(jié)果生成了每一條染色體等位基因頻率分布pdf文件

9.png
x_11.png
x_12.png
x_13.png
x_14.png
  • SHOREmap annotate：注釋候選區(qū)域
SHOREmap annotate --chrsizes chromsize.txt --folder ../analysis/ --snp ../foreground/1_converted_variant.txt --chrom Chr3 --start 1 --end 4000000 --genome ../genome/TAIR10_chr_all.fas -gff../annotation/TAIR10_GFF3_genes.gff
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容