作者：rapunzel0103

鏈接：http://www.reibang.com/p/53362fe881cd

最近跑通了一遍GWAS分析，全程在linux操作，雖然具體還有好多需要微調(diào)的地方繁仁，先把代碼整理分享出來mark一下

前期準(zhǔn)備

1.理論知識(shí)

強(qiáng)烈推薦百邁客云課堂課程GWAS生物信息培訓(xùn)課程

或者可以看看我的gwas相關(guān)文章

GWAS基本分析內(nèi)容（課程學(xué)習(xí)筆記）

常用GWAS統(tǒng)計(jì)方法和模型簡(jiǎn)介（課程學(xué)習(xí)筆記）

臨床生物信息學(xué)中的GWAS分析（內(nèi)附擴(kuò)展閱讀）

精細(xì)定位——降低 GWAS的復(fù)雜度（文獻(xiàn)研讀）

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

如果你沒有自己的數(shù)據(jù)又想做gwas分析的話，可以選擇3000水稻基因組的http://snp-seek.irri.org/數(shù)據(jù)庫直接下載，vcf、表型數(shù)據(jù)甚至plink bed/bim/fam文件直接下載级遭，gwas結(jié)果也做到了可視化

另外推薦華農(nóng)謝為博團(tuán)隊(duì)開發(fā)的這個(gè)網(wǎng)站http://ricevarmap.ncpgr.cn/v2/非常好用，提供數(shù)據(jù)下載和gwas可視化結(jié)果）

基因型數(shù)據(jù)可以根據(jù)bioproject accession編號(hào)從NCBI上下載：

表型數(shù)據(jù)直接從網(wǎng)站上以excel或csv格式導(dǎo)出：

533份樣品下載和比對(duì)挺耗時(shí)和占內(nèi)存的渺尘，建議保留bam文件挫鸽，其他的包括fastq（由bam能轉(zhuǎn)成fq）、中間文件都可以刪掉鸥跟，數(shù)據(jù)過濾質(zhì)控很快丢郊，整個(gè)項(xiàng)目做下來大概耗時(shí)一個(gè)多月吧

第一步 SNP calling

需要安裝的軟件：BWA和GATK/samtools

一盔沫、BWA比對(duì)

1.構(gòu)建index

bwa index -a is ref.fa ?或bwa index -a bwtsw ref.fa (>2G)

samtools faidx ref.fa

java -jar $picard/CreatSequenceDictionary.jar R=ref.fa O=ref.dict

2.每個(gè)樣分別比對(duì)到參考基因組

bwa mem? -t 5? -M -R "@RG\tID:A\tSM:A" ref.fa? A1_1.fq A1_2.fq > A1.sam ?&

bwa mem? -t 5? -M -R "@RG\tID:A\tSM:A" ref.fa? A2_1.fq A2_2.fq > A2.sam ? & ? ? 以此類推.......（-M 將shorter split hits標(biāo)記為次優(yōu)，可以兼容Picard．-R 每個(gè)標(biāo)記號(hào)需不同枫匾，方便后面合并）

3.SortSam

java -jar $picard/SortSam.jar I=A.sam O=A1.sort.bam SO=coordinate

4.MarkDuplicates

java -jar $picard/MarkDuplicates.jar I=A1.sort.bam O=A1.Mark.bam M=A1.metrics?

二架诞、SNP檢測(cè)

1.RealignerTargetCreator

java -jar $GATK -R $ref -T RealignerTargetCreator -I A1.Mark.bam -o A1.realign.interval_list

2.IndelRealigner

java -jar $GATK -R $ref -T IndelRealigner -I A1.Mark.bam? -o A1.realn.bam -targetIntervals A1.realign.interval_list

3.HaplotypeCaller

?java -jar $GATK -T HaplotypeCaller -R $ref -ERC GVCF -I A1.realn.bam? ? ?? --variant_index_type LINEAR --variant_index_parameter 128000 ? ? -o A1.gvcf

4.CombineGVCFs

?java -jar $GATK -T CombineGVCFs -R $ref --disable_auto_index_creation_and_locking_when_reading_rods --variant A1.gvcf?A2.gvcf?A3.gvcf??.... ? -o combine.gvcf ?(這一步需要把每個(gè)樣品的gvcf合并)

5.GenotypeGVCFs

java?? -jar $GATK -T GenotypeGVCFs -nt 4 -R $ref --disable_auto_index_creation_and_locking_when_reading_rods? ?? -o test_final.vcf --variant combine.gvcf

第二步 基因型填補(bǔ)

需要安裝的軟件：Tassel/beagle等

1. tassel (-LDKNNilmputationPlugin參數(shù)有誤，沒有跑成功, 有朋友指出LDKNNi后面是大寫的i干茉，不是小寫的L,以后再試試)

perl /home/user/soft/tassel_v5/run_pipeline.pl -Xms512m -Xmx5g -importGuess ?test_final.vcf -LDKNNiImputationPlugin -highLDSSites 50 -knnTaxa 10 -maxLDDistance 100000000 -endPlugin -export test.imputed.vcf -exportType VCF

也可以java -jar sTASSEL.jar 在窗口操作 記得給服務(wù)器接顯示屏

2.beagle

java -jar beagle.08Jun17.d8b.jar gt=test_final.vcf?out=test.imputed.vcf

第三步 數(shù)據(jù)篩選及格式轉(zhuǎn)換

需要安裝的軟件：plink等

1.按MAF>0.05和缺失率<0.1過濾

/home/user/soft/plink --vcf test.imputed.vcf--maf 0.05 --geno 0.1--recode vcf-iid --out test.filter --allow-extra-chr （非數(shù)字染色體號(hào)ChrUn/Sy用此參數(shù), 建議盡量把染色體號(hào)轉(zhuǎn)成數(shù)字谴忧，另外需要對(duì)vcf中的標(biāo)記ID進(jìn)行編號(hào)）

2.對(duì)標(biāo)記進(jìn)行LD篩選

/home/user/soft/plink --vcf test.filter.vcf--indep-pairwise 50 10 0.2--out test.filterLD --allow-extra-chr （.in文件里是入選的標(biāo)記id）

3.提取篩選結(jié)果

/home/user/soft/plink --vcf test.filterLD.vcf --make-bed--extracttest.filter.in --out? test.filter.prune.in

4.轉(zhuǎn)換成structure/admixture格式

/soft/plink --bfile test.filter.prune.in--recodestructure--out test.filter.prune.in ?#生成. recode.strct_in為structure輸入格式

/soft/plink --bfile test.filter.prune.in--recode12--out test.filter.prune.in ?#生成.ped為admixture輸入格式

第四步 群體結(jié)構(gòu)

需要安裝的軟件：structure/admixture等

這里我選了比較簡(jiǎn)單admixture來做 k值范圍1到13

/soft/admixture_linux-1.3.0/admixture --cv test.filter.ped 1 >>log.txt

/soft/admixture_linux-1.3.0/admixture --cv test.filter.ped 2 >>log.txt

.......

/soft/admixture_linux-1.3.0/admixture --cv test.filter.ped 13 >>log.txt

wait

grep "CV error" log.txt >k_1to13

取CV error最小時(shí)的k值=10, ?其中test.filter.prune.in.10.Q結(jié)果文件作為關(guān)聯(lián)分析的輸入源文件（去掉最后一列 添加表頭和ID）

第五步 親緣關(guān)系/PCA(選做)

需要安裝的軟件: tassel等（ PCA分析可以用R包，已經(jīng)做了群體結(jié)構(gòu)這里就沒做PCA分析）

perl?/tassel_v5/run_pipeline.pl -importGuess test_impute.vcf -KinshipPlugin ?-method Centered_IBS -endPlugin -export test_kinship.txt -exportType SqrMatrix

第六步 關(guān)聯(lián)分析

需要安裝的軟件: tassel/GAPIT/FaSt-LMM等（ GAPIT很強(qiáng)大等脂，要裝很多R包俏蛮，自動(dòng)做圖可視化撑蚌。FaSt-LMM以后打算嘗試一下）

輸入文件的格式需要手動(dòng)修改一下上遥，也比較簡(jiǎn)單 (如圖)

test.best_k10.txt

test.trait.txt

test_kinship.txt

1.vcf轉(zhuǎn)hapmap格式

perl /tassel_v5/run_pipeline.pl -fork1 -vcf test.imputed.vcf? -export test -exportType Hapmap -runfork1

2.SNP位點(diǎn)排序

perl /tassel_v5/run_pipeline.pl -SortGenotypeFilePlugin -inputFile test.hmp.txt? -outputFile test_sort -fileType Hapmap

得到test.hmp.txt

3. GLM模型

perl /tassel_v5/run_pipeline.pl -fork1 -h test_sort.hmp.txt -fork2 -r test.trait.txt -fork3 -q test.best_k10.txt -excludeLastTrait -combine4 -input1 -input2 -input3 -intersect -glm -export test_glm -runfork1 -runfork2 -runfork3

4.MLM模型

perl /tassel_v5/run_pipeline.pl -fork1 -h test_sort.hmp.txt -fork2 -r test.trait.txt -fork3 -q test.best_k10.txt -excludeLastTrait -fork4 -k test_kinship.txt -combine5? -input1 -input2 -input3 -intersect -combine6 -input5 -input4 -mlm -mlmVarCompEst P3D -mlmCompressionLevel None -export test_mlm -runfork1 -runfork2 -runfork3

最后得到關(guān)聯(lián)分析的結(jié)果文件，喜大普奔争涌！

當(dāng)然GWAS分析需要根據(jù)實(shí)際項(xiàng)目材料的需要粉楚，靈活地選擇分析方法，理解統(tǒng)計(jì)學(xué)亮垫、群體遺傳學(xué)等原理模软，認(rèn)識(shí)GWAS的特點(diǎn)和局限性還是很有必要的。這里只是簡(jiǎn)單地在linux上跑通了一遍常用的流程饮潦，有很多R包MVP燃异，GAPIT等等可以做并且一鍵出圖如CMplot～有很多不懂的地方需要多學(xué)習(xí)，畢竟大牛都是自己手動(dòng)寫腳本來分析继蜡，不怎么用軟件的回俐。。/(ㄒoㄒ)/

其他參考：http://blog.sina.com.cn/s/blog_83f77c940102w3d2.html

https://wenku.baidu.com/view/efdb4115a2161479171128e9.html

作者：rapunzel0103

鏈接：http://www.reibang.com/p/53362fe881cd

來源：簡(jiǎn)書

著作權(quán)歸作者所有稀并。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)仅颇，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處。