偶然發(fā)現(xiàn)有個(gè)R包可以做相關(guān)的QTL 分析即寡。拿來試試畦木。
首先 安裝,參考官方文檔https://github.com/bmansfeld/QTLseqr.
第一步萍桌，如果沒有安裝devtools的話古话，install.packages("devtools")牺氨，有就直接調(diào)用，如果這個(gè)不會草娜，就不用往下看了挑胸。
第二步，devtools::install_github("bmansfeld/QTLseqr")
從github加載QTLseqr宰闰。
拿他的論文實(shí)例來演練一次茬贵。
library("QTLseqr")
看一下輸入文件，文檔中說不僅可以識別GATK中VariantsToTable 功能下的table格式的文件移袍，也能識別包含每個(gè)批量的等位基因讀取深度的分隔文件解藻。調(diào)用的功能是importFromGATK and mportFromTable。
我們先從GitHub下載一個(gè)Yang2013data 葡盗。
命令如下：devtools:: install_github("bmansfeld/Yang2013data")
library("Yang2013data")?
導(dǎo)入數(shù)據(jù)
rawData <-system.file(
"extdata",
"Yang_et_al_2013.table",
package ="Yang2013data",
mustWork =TRUE)
如果自己有數(shù)據(jù)table格式：
rawdata<-(file_to_path/my_table_rawData <-system.file(
"extdata",
"Yang_et_al_2013.table",
package ="Yang2013data",
mustWork =TRUE)
如果自己有數(shù)據(jù)
rawdata<-(file_to_path/my_table)
然后我們分別給高低池命名螟左。寫染色體編號。
HighBulk <-"SRR834931"
LowBulk <-"SRR834927"
Chroms <-paste0(rep("Chr",12),1: 12
用GATK call snp觅够。table文件用線面這個(gè)命令轉(zhuǎn)換胶背。
java -jar GenomeAnalysisTK.jar \
-T VariantsToTable \
-R ${REF} \
-V ${NAME} \
-F CHROM -F POS -F REF -F ALT \
-GF AD -GF DP -GF GQ -GF PL \
-o ${NAME}.table
看看ImportFromGATK 功能
df <-
importFromGATK(
file =rawData,
highBulk =HighBulk,
lowBulk =LowBulk,
chromList =Chroms
)
會計(jì)算兩個(gè)池的等位基因頻率，snp-index喘先，還有Δindex钳吟。
導(dǎo)入一下分隔文件【秸可以是csv, tsv或其他格式红且。格式表頭如下。
如果有現(xiàn)成的table文件涤姊，也可以importFromTable
df <-importFromTable(file ="Yang2013.csv",
highBulk =HighBulk,lowBulk =LowBulk,
chromList =Chroms)
矩陣如下：## CHROM POS REF ALT AD_REF.LOW AD_ALT.LOW DP.LOW GQ.LOW PL.LOW
## 1 Chr1 31071 A G 34 36 70 99 897,0,855
## 2 Chr1 31478 C T 34 52 86 99 1363,0,844
## 3 Chr1 33667 A G 20 48 68 99 1331,0,438
## 4 Chr1 34057 C T 38 40 78 99 1059,0,996
## 5 Chr1 35239 A C 25 36 61 99 987,0,630
## 6 Chr1 38389 T C 36 42 78 99 1066,0,906
## SNPindex.LOW AD_REF.HIGH AD_ALT.HIGH DP.HIGH GQ.HIGH PL.HIGH
## 1 0.5142857 26 22 48 99 522,0,698
## 2 0.6046512 40 34 74 99 848,0,1099
## 3 0.7058824 24 29 53 99 765,0,599
## 4 0.5128205 29 26 55 99 673,0,772
## 5 0.5901639 40 60 100 99 1632,0,1015
## 6 0.5384615 42 40 82 99 984,0,1105
## SNPindex.HIGH REF_FRQ deltaSNP
## 1 0.4583333 0.5084746 -0.055952381
## 2 0.4594595 0.4625000 -0.145191703
## 3 0.5471698 0.3636364 -0.158712542
## 4 0.4727273 0.5037594 -0.040093240
## 5 0.6000000 0.4037267 0.009836066
## 6 0.4878049 0.4875000 -0.050656660
先來看看測序深度的直方圖
再來看看參考基因組的頻率分布
再來看一下每個(gè)混池的snp-index暇番，對于F2群體，大部分位點(diǎn)再0.5.（非連鎖位點(diǎn)）
過濾SNP位點(diǎn)
df_filt <-
filterSNPs(
SNPset =df,
refAlleleFreq =0.20,
minTotalDepth =100,
maxTotalDepth =400,
depthDifference =100,
minSampleDepth =40,
minGQ =99,
verbose =TRUE
)
QTLseq analysis
開始分析
df_filt <-runQTLseqAnalysis(df_filt,
windowSize =1e6,
popStruc ="F2",
bulkSize =c(385,430),
replications =10000,
intervals =c(95,99)
)
這里好像只適用于F2與RIL思喊，這個(gè)世界總是對BIL充滿了惡意壁酬。沒關(guān)系，BIL親測也適用，
G分析厨喂。
df_filt <-runGprimeAnalysis(df_filt,
windowSize =1e6,
outlierFilter ="deltaSNP",
filterThreshold =0.1)
現(xiàn)在開始QTL分析
現(xiàn)在我們對過濾后的數(shù)據(jù)感到滿意，看起來G ’接近于對數(shù)正態(tài)分布庄呈，我們
可以繪制一些全基因組數(shù)據(jù)并最終嘗試識別QTL蜕煌。
p1 <-plotQTLStats(SNPset =df_filt,var ="nSNPs")
從p1圖中，我們可以得到全基因組snp密度诬留。
p2 <-plotQTLStats(SNPset =df_filt,var ="deltaSNP",plotIntervals =TRUE)
超過置信區(qū)間的有4個(gè)斜纪。其中紅色p<0.05,綠色P<0.01.
我們再來看一下G'是否是顯著的，并且FDR (q)<0.01
p3 <-plotQTLStats(SNPset =df_filt,var ="Gprime",plotThreshold =TRUE,q =0.01)
QTLplots <-plotQTLStats(
SNPset =df_filt,var ="negLog10Pval",
plotThreshold =TRUE,
q =0.01,
subset =c("Chr1","Chr8")
)
提取顯著性區(qū)間的信息
QTL <-getSigRegions(SNPset =df_filt,alpha =0.01)
提取出csv格式的數(shù)據(jù)
results <-getQTLTable(SNPset =df_filt,method ="Gprime",alpha =0.01,export =FALSE)
