使用SHOREmap做mapping-by-sequencing

簡介

SHOREmap可以用來分析傳統(tǒng)作圖群體（自然系natural strains和分化系,diverged accession雜交全蝶，或outcrossing)或近等作圖群體（isogenic mapping population, 誘變后代與未誘變親本進行雜交，即會交,backcrossing)所產(chǎn)生的重測序數(shù)據(jù)。根據(jù)作圖群體構(gòu)建方式的不同敏释，SHOREmap的outcross或backcross采用不同基于滑窗(sliding)方式對等位基因頻率進行分析捞烟。

SHOREmap的backcross和outcross都需要從突變重組庫中獲得的一致的堿基識別信息

安裝

前置安裝

SHOREmap需要DISLIN科學(xué)庫進行數(shù)據(jù)可視化

但是在安裝DISLIN之前還需要保證存在/usr/lib/libXm.so*和/usr/lib/libXm.so*嗽仪，這兩者的安全需要root權(quán)限，所以要么聯(lián)系管理員俩檬，要么想辦法繞開(這個辦法，我還沒有想到).

sudo apt-get update
sudo apt-get install libmotif4
sudo apt-get install libxt-dev

開始安裝doslin庫

cd /path/to/src
# 下載
wget ftp://ftp.gwdg.de/pub/grafik/dislin/linux/i586_64/dislin-11.0.linux.i586_64.tar.gz
# 解壓縮
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 ./example/dislin_d.h $DISLIN
# 刪除安裝文件(可選)
rm -rf dislin-11.0

安裝SHOREmap v3.x

我這次安裝的是當(dāng)前最新的3.4版本碾盟，其他版本估計換湯不換藥棚辽。

cd $HOME/biosoft
wget http://bioinfo.mpipz.mpg.de/shoremap/SHOREmap_v3.4.tar.gz
# 替換SHOREmap下的dislin的一些文件
tar -zxvf SHOREmap_v3,4
rm dislin/*dislin_d.*
cp $DISLIN/*dislin_d.* dislin
# 編輯/etc/profile或.bashrc
vi .bashrc
export LD_LIBRARY_PATH=$HOME/src/SHOREmap_v3.4/dislin
# 退出保存.bashrc: Esc+:wq
source .bashrc
# 到之前安裝的文件夾下
cd & cd src/SHOREmap_v3.4
(可選，如果沒有g(shù)++)sudo apt-get install build-essential
make
# 將編譯文件拷貝到習(xí)慣的文件夾中巷疼，然后添加執(zhí)行路徑
cp SHOREmap ../../biosoft/SHOREmap_v3.4
echo "export $HOME/bisoft/SHOREmap_v3.4" >> ~/.bashrc

最后晚胡，可以重新啟動一下bash驗證

官方網(wǎng)站提供的兩個常見問題的解答

Note 1: if the compilation complains like "/usr/bin/ld: cannot find -lXt" (or "/usr/bin/ld: cannot find -ldislin_d"), please open the makefile with the command

vi makefile

Press keys 'esc' and 'i' on the keyboard to edit makefile; move the cursor with arrow keys to the position before -lXt, and edit -L/path/to/libXt.so/; if '-ldislin_d' is not found, edit -L/path/to/dislin_d/ before -ldislin_d). After that, press keys 'esc', type in :wq, and press enter to save editing and quit vi. ('-L' tells the linker where to find the library given by -l)

Note 2: if '/usr/lib/ld: warning: libXm.so.3, needed by ./dislin/libdislin_d.so, not found (try using -rpath or -rpath-link)' occurs, and you have installed libmotif4, do the following:

cp /usr/lib/libXm.so.4 /usr/lib/libXm.so.3

We can make SHOREmap avaiable for general use by inserting the following command into /etc/profile

export PATH=$PATH:/path/to/SHOREmap_v3.x

and

source /etc/profile

總體流程

OUTCROSS

outcross的基本步驟 描述
SHOREmap extract 提取與SNP突變相關(guān)的重測序的一致的識別
SHOREmap create 根據(jù)背景/親本系的重測序質(zhì)量創(chuàng)建SNP標(biāo)記列表
SHOREmap outcross 進行等位基因頻率分析并定義mapping interval（也就是找到突變所在的大致區(qū)域)
SHOREmap annotate 對mapping interval中的突變基因效應(yīng)進行注釋

BACKCROSS

backcross的基本步驟 描述
SHOREmap extract 提取與SNP突變相關(guān)的重測序的一致的識別
SHOREmap backcross 進行等位基因分析
SHOREmap annotate 對mapping interval中的突變基因效應(yīng)進行注

具體步驟

下載數(shù)據(jù)

只安裝軟件，卻沒有數(shù)據(jù)嚼沿，我們也只能干瞪眼估盘。

oucross分析所需數(shù)據(jù)

# OCF2 
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.fg.reads1.fq.gz &
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.fg.reads2.fq.gz &
# Ler 
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.bg.reads1.fq.gz &
wget -4 -qh ttp://bioinfo.mpipz.mpg.de/shoremap/data/software/OC.bg.reads2.fq.gz &

backcross分析所需數(shù)據(jù)

# BCF2
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads1.fq.gz &
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.fg.reads2.fq.gz &
# mir159a (Col)
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads1.fq.gz
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/BC.bg.reads2.fq.gz

其他數(shù)據(jù)

除了最基本的測序數(shù)據(jù)外，我們可能還需要參考基因組骡尽，已有的注釋數(shù)據(jù)等

# 參考基因組
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/TAIR10_chr_all.fas &
# 基因注釋
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/TAIR10_GFF3_genes.gff &
# SHORE操作的結(jié)果數(shù)據(jù)
wget -4 -q http://bioinfo.mpipz.mpg.de/shoremap/data/software/scoring_matrix_het.txt &

重測序

首先使用bwa,bowtie2等read比對工具將得到的數(shù)據(jù)比對到參考基因組上遣妥。
假設(shè)你當(dāng)前處在MBS文件夾下，該文件下有如下文件

# 混池測序結(jié)果攀细，雙端
BC.fg.reads1.fq.gz
BC.fg.reads2.fq.gz
# 背景信息測序結(jié)果
BC.bg.reads1.fq.gz
BC.bg.reads2.fa.gz
# 擬南芥參考基因組
TAIR10_chr_all.fas
# 擬南芥注釋信息
TAIR10_GFF3_genes.gff

以下操作都是基于上述文件進行箫踩。

第一步: 序列比對，產(chǎn)生SAM文件

mkdir index
# 創(chuàng)建比對所需索引
bowtie2-build TAIR10_chr_all.fas index/TAIR10
# 序列比對
bowtie2 -x index/TAIR10 -1 BC.fg.reads1.fq.gz -2 BC.fg.reads2.fq.gz -S FG.sam
bowtie2 -x index/TAIR10 -1 BC.bg.reads1.fq.gz -2 BC.bg.reads2.fq.gz -S BG.sam

第二步： SAMtools預(yù)測突變位點

為了加快運算速度谭贪，可以先轉(zhuǎn)換格式境钟，并排序

samtools view -b -o BG.bam BG.sam
samtools view -b -o FG.bam FG.sam
samtools sort -o BG.sorted.bam BG.bam 
samtools sort -o FG.sorted.bam FG.bam 
samtools index BG.sorted.bam
samtools index FG.sorted.bam

consensus-calling program 尋找可能的變異位點
由于官方說明中使用samtools版本為0.1.19，先要解釋一下參數(shù) mpileup -uD -f俭识，

# -f:faidx indexed reference sequence file 前后版本一致
# -u:generate uncompress BCF output 前后版本一致
# -D:output per-sample DP in BCF (require -g/-u).與輸出格式有關(guān)慨削，目前改為-t

因此，對于1.4版本的samtools套媚，相應(yīng)的參數(shù)為

mpileup -u -t DP -f

官方說明的bcftools也是0.1.19缚态，參數(shù)為bcftools view -vcg 舊版本的view在當(dāng)前的版本用于過濾，功能被call替代

# -v Output variant sites only (force -c)
# -c屬于Call variants using Bayesian inference. This option automatically invokes option -e.When -v is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0],目前被m取代
# -g Call per-sample genotypes at variant sites (force -c)堤瘤，這個沒有找到合適的替代參數(shù)

綜上玫芦，推薦使用如下命令行

samtools mpileup -u -t DP -f ../../../index/TAIR10_chr_all.fa ../../align/bwa/default/BG.sorted.bam | bcftools call -vm -Ov > BG.vcf &
samtools mpileup -u -t DP -f ../../../index/TAIR10_chr_all.fa ../../align/bwa/default/FG.sorted.bam | bcftools call -vm -Ov > FG.vcf &

額外步驟：VCF格式轉(zhuǎn)換

由于vcftools工具版本，所以最后的文件版本是4.2本辐，而SHOREmap要求4.1桥帆。通過biostar找到高人寫的降級工具（其實就是把一些字符替換一下，但是不了解vcf不同版本的差異話慎皱，是不知道怎么寫）
把下面的代碼存為vcf_dowgrade.sh

# If you are trying to view VCF 4.2 files in IGV - you may run into issues. This function might help you.
# This script will:
# 1. Rename the file as version 4.1
# 2. Replace parentheses in the INFO lines (IGV doesn't like these!)

function vcf_downgrade() {
  outfile=${1/.bcf/}
  outfile=${outfile/.gz/}
  outfile=${outfile/.vcf/}
  bcftools view --max-alleles 2 -O v $1 | \
  sed "s/##fileformat=VCFv4.2/##fileformat=VCFv4.1/" | \
  sed "s/(//" | \
  sed "s/)//" | \
  sed "s/,Version=\"3\">/>/" | \
  bcftools view -O z > ${outfile}.dg.vcf.gz
  tabix ${outfile}.dg.vcf.gz
}

其實對于單個文件而言环葵，可以直接用以下命令

infile=BG.vcf
outfile=BG.vcf
bcftools  view --max-alleles 2 -O v ${infile} | \
sed "s/##fileformat=VCFv4.2/##fileformat=VCFv4.1/" | \
  sed "s/(//" | \
  sed "s/)//" | \
  sed "s/,Version=\"3\">/>/" | \
  bcftools view -O z > ${outfile}.dg.vcf.gz

使用SHOREmap尋找突變所在區(qū)

第一步：需要把bcf文件通過SHOREmap convert轉(zhuǎn)換成SHOREmap能認(rèn)識的格式

SHOREmap convert --marker samtools.vcf --folder path/to/folder -runid int

會生成三個文件3_converted_consen.txt, 3_converted_variant.txt and 3_converted_reference.txt.

第二步：提取候選分子標(biāo)記的consensus information(mapping pool)

SHOREmap extract --chrsizes chromsize.txt --folder ../SHOREmap_analysis --marker 11_converted_variant.txt --consen 11_converted_consen.txt -verbose

第三步：
然后使用SHOREmap backcross分析。SHOREmap backcross可用來分析回交作圖群體所得到重組后代混池數(shù)據(jù)宝冕。相對于傳統(tǒng)作圖群體，只有誘變劑產(chǎn)生的突變會分離邓萨，也只有這些才會用于突變定位地梨。

SHOREmap backcross會嘗試過濾出所有參考基因組和測序池之間不同部分用于找到突變點特異部分菊卷。為了保證不是自然變異或者是測序錯誤，測序池選擇的部分要多次出現(xiàn)在親本或背景中宝剖。然后根據(jù)前景和/或背景的(識別堿基,base calls,質(zhì)量/覆蓋率/等位基因)信息洁闰，確定是否把保留的SNP位點作為分子標(biāo)記。在正確的篩選后（擬南芥大概有上百個標(biāo)記）万细，SHOREmap backcross就能在分析marker的AF后識別大致的峰扑眉。進一步對變異注釋后，就能找到目標(biāo)性狀的候選基因了赖钞。

SHOREmap backcross所需的輸入文件如下：

  • 染色體大小文件腰素，--chrsizes。分為兩行雪营，一行是染色體位置弓千，一行是染色體大小。scaffold同理
  • 候選marker文件献起。也就是使用SHOREmap convert通過vcf生成的converted_variant.txt洋访，每一列的含義如下。
    Column Description
    1 Project name
    2 Identity of chromosome
    3 Position of the SNP-marker
    4 Reference base
    5 Alternative base (or mutant base)
    6 Quality of the alternative base (ranging from 0 to 40)
    7 Number of reads supporting the predicted base
    8 Ratio of reads supporting the predicted base to total coverage
SHOREmap backcross --chrsizes chromsize.txt --marker ../convert/11_converted_variant.txt --consen extracted_consensus_0.txt --folder ../BC_analysis -plot-bc --marker-score 40 --marker-freq 0.0 --min-converage 10 --max-coverage 80 -bg ../convert/12_converted_variant.txt  --bg-cov 1 --bg-freq 0.0 --bg-score 1 -non-EMS --cluster 1 --marker-hit 1 -verbose

第四步：對結(jié)果進行注釋

SHOREmap annotate --chrsizes chromsize.txt --folder ../BC_analysis/ann --snp ../convert/11_converted_variant.txt --chrom 2 --start 1 --end 4000000 --genome ../../TAIR10_chr_all.fas --gff ../../TAIR10_GFF3_genes.gff

注意點：

Alignment：不同的比對軟件和參數(shù)設(shè)計會對結(jié)果有多大影響
SNP calling: samtools+vcftools 或者是GATK產(chǎn)生的vcf文件所包含的結(jié)果對結(jié)果有什么影響
SHOREmap： 如何調(diào)整閾值谴餐，如何辨認(rèn)結(jié)果姻政。

目標(biāo)：

比對基本上是高通量的數(shù)據(jù)分析的必經(jīng)之路，必須要對比對有足夠的理解岂嗓。
vcf格式描述了

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末汁展，一起剝皮案震驚了整個濱河市，隨后出現(xiàn)的幾起案子摄闸，更是在濱河造成了極大的恐慌善镰，老刑警劉巖，帶你破解...
    沈念sama閱讀 222,807評論 6贊 518
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件年枕，死亡現(xiàn)場離奇詭異炫欺，居然都是意外死亡，警方通過查閱死者的電腦和手機熏兄，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,284評論 3贊 399
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進店門品洛，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人摩桶，你說我怎么就攤上這事桥状。” “怎么了硝清？”我有些...
    開封第一講書人閱讀 169,589評論 0贊 363
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵辅斟，是天一觀的道長。 經(jīng)常有香客問我芦拿，道長士飒，這世上最難降的妖魔是什么查邢？ 我笑而不...
    開封第一講書人閱讀 60,188評論 1贊 300
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮酵幕，結(jié)果婚禮上扰藕，老公的妹妹穿的比我還像新娘。我一直安慰自己芳撒，他們只是感情好邓深，可當(dāng)我...
    茶點故事閱讀 69,185評論 6贊 398
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布。 她就那樣靜靜地躺著笔刹，像睡著了一般芥备。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上徘熔，一...
    開封第一講書人閱讀 52,785評論 1贊 314
  • 城市分裂傳說
    那天门躯，我揣著相機與錄音，去河邊找鬼酷师。 笑死讶凉，一個胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的山孔。 我是一名探鬼主播懂讯，決...
    沈念sama閱讀 41,220評論 3贊 423
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼台颠！你這毒婦竟也來了褐望？” 一聲冷哼從身側(cè)響起，我...
    開封第一講書人閱讀 40,167評論 0贊 277
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤串前，失蹤者是張志新（化名）和其女友劉穎瘫里，沒想到半個月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體荡碾，經(jīng)...
    沈念sama閱讀 46,698評論 1贊 320
  • ?護林員之死
    正文 獨居荒郊野嶺守林人離奇死亡谨读，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 38,767評論 3贊 343
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了坛吁。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片劳殖。...
    茶點故事閱讀 40,912評論 1贊 353
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖拨脉，靈堂內(nèi)的尸體忽然破棺而出哆姻，到底是詐尸還是另有隱情，我是刑警寧澤玫膀，帶...
    沈念sama閱讀 36,572評論 5贊 351
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布矛缨，位于F島的核電站，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏劳景。R本人自食惡果不足惜誉简，卻給世界環(huán)境...
    茶點故事閱讀 42,254評論 3贊 336
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望盟广。 院中可真熱鬧，春花似錦瓮钥、人聲如沸筋量。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,746評論 0贊 25
  • 一樁弒父案碉熄，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽桨武。三九已至，卻和暖如春锈津，著一層夾襖步出監(jiān)牢的瞬間呀酸，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,859評論 1贊 274
  • 情欲美人皮
    我被黑心中介騙來泰國打工琼梆， 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留性誉，地道東北人。 一個月前我還...
    沈念sama閱讀 49,359評論 3贊 379
  • 代替公主和親
    正文 我出身青樓茎杂，卻偏偏與公主長得像错览，于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子煌往，可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 45,922評論 2贊 361

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

  • Spring Cloud為開發(fā)人員提供了快速構(gòu)建分布式系統(tǒng)中一些常見模式的工具（例如配置管理倾哺，服務(wù)發(fā)現(xiàn)，斷路器刽脖，智...
    卡卡羅2017閱讀 134,716評論 18贊 139
  • 吸管只能用來喝水?孩子們的腦洞那么大羞海，怎么可能放過任何好玩的東西，在游戲中發(fā)揮創(chuàng)意曲管，在故事里實現(xiàn)夢想却邓，在生活中感悟...
    設(shè)計獅媽咪閱讀 1,481評論 0贊 0
  • 如果你在一個企業(yè)學(xué)不到自己想要的知識衙耕，工作經(jīng)驗半年昧穿，你會跳槽嗎？ 被問到這個問題橙喘，感覺回到了當(dāng)初的自己时鸵。那時候，我...
    左手清香閱讀 174評論 0贊 0
  • 夢到經(jīng)歷了很多很多的事情，最后父親用電瓶車載我回家饰潜，醒來后鼻子一酸初坠，真是對不起自己的家人，對不起自己逝去的歲月
    268fb13a2916閱讀 264評論 0贊 0
  • “嘭”的一聲半沽，橙子先生的車門被關(guān)上，接著就看到一個清瘦的身影漸漸走遠(yuǎn)了吴菠，沒有回頭者填。 誰也不知道安妮心中此時在想些...
    等風(fēng)來cll閱讀 325評論 0贊 2