直行依序為:1.各株外表型data、2~47:每一株的genotype data
最後把檔案存成csv檔(我直接存在D槽)
(二)資料輸入
下載package"qtl"
>library(qtl)
> test01<-read.cross("csv",dir="d:",file="RqtlLA1543.csv")
# dir是指資料存放處可以先用?read.cross看其用法!!
--Read the following data:
46 individuals
122 markers
1 phenotypes
--Estimating genetic map
--Cross type: f2
Warning message:
In summary.cross(cross) :
Some chromosomes > 1000 cM in length; there may be a problem with the genetic map.
(Perhaps it is in basepairs?)
> test01
This is an object of class "cross".
It is too complex to print, so we provide just this summary.
F2 intercross
No. individuals: 46
No. phenotypes: 1
Percent phenotyped: 100
No. chromosomes: 12
Autosomes: 1 2 3 4 5 6 7 8 9 10 11 12
Total markers: 122
No. markers: 15 13 12 10 11 7 11 8 7 9 11 8
Percent genotyped: 96.9
Genotypes (%): AA:22.3 AB:50.2 BB:27.4 not BB:0 not AA:0
> formLinkageGroups(test01)
#重新分配連鎖群，可推估還需要添加那些marker。
origchr LG #origchr:在data中已先assigned的連鎖群，LG重新估算後的。
SSR98 1 1
M02 1 1
SSR316 1 1
SSR134 1 1
M4 1 1
M05 1 1
M06 1 1
SSR595 1 1
M07 1 1
TM10 1 1
TM110 1 1
T1-5 1 1
T1-9 1 1
T1-7 1 1
M09 2 2
SSR32 2 2
SSR331 2 2
M10 2 2
T2-11 2 2
T2-1 2 2
T2-2 2 2
T2-3 2 2
T2-8 2 2
T2-10 2 2
T2-12 2 2
Ovtcaps 2 2
OvtEx2 2 2
M11 3 3
M12 3 3
M15 3 3
M16 3 3
M17 3 3
TM40 3 3
T3-1 3 3
T3-3 3 3
T3-2 3 3
T3-4 3 3
T3-5 3 3
T3-6 3 3
M19 4 7
M20 4 7
M21 4 14
T4-1 4 14
T4-2 4 7
T4-3 4 7
T4-4 4 7
T4-5 4 7
T4-6 4 7
T4-8 4 7
> test02<-est.rf(test01)
> plotRF(test02)
#計算LOD值，即可知道QTL位置應該位於何處。但scanone僅適合單一QTL的找尋。
M10 2 8.426426 4.22771445
T2-10 2 13.098849 5.27987801
OvtEx2 2 20.198443 5.36025217
Ovtcaps 2 22.551423 5.33027185
SSR331 2 30.689005 4.65098199
SSR32 2 36.517105 3.88908267
T2-2 2 46.713341 3.10885368
M65 12 0.000000 2.99672232
M66 12 5.743001 3.55032910
T12-1 12 21.974863 3.43095251
M67 12 28.863282 1.69070880
SSR345 12 44.236000 0.76826234
Tg50 12 55.627526 1.06820054
T12-3 12 73.950488 0.41372156
M68 12 90.821113 0.20482765
> test02<-est.rf(test01)
> plotRF(test02)
#計算LOD值，即可知道QTL位置應該位於何處。但scanone僅適合單一QTL的找尋。
> scanone(test01)
chr pos lod
SSR98 1 0.000000 0.45531293
TM110 1 0.000000 0.45531293
SSR316 1 11.733378 0.31931487
SSR134 1 23.906352 0.27656684
M02 1 25.074505 0.17178144
T1-5 1 32.863114 1.46767194
T1-9 1 36.679146 1.41270063
M4 1 49.130293 1.76414171
T1-7 1 56.364825 1.01114976
M05 1 67.360089 0.23521435
M06 1 73.011499 0.11050005
SSR595 1 83.695745 0.09456710
M07 1 97.952843 0.11867296
TM10 1 101.444432 0.16816227
T2-1 2 0.000000 3.07644668M10 2 8.426426 4.22771445
T2-10 2 13.098849 5.27987801
OvtEx2 2 20.198443 5.36025217
Ovtcaps 2 22.551423 5.33027185
SSR331 2 30.689005 4.65098199
SSR32 2 36.517105 3.88908267
T2-2 2 46.713341 3.10885368
T2-3 2 55.830503 2.18501681
M09 2 60.706826 2.30268401
T2-12 2 76.012547 0.34917303
T2-11 2 91.992752 0.04704214
M65 12 0.000000 2.99672232
M66 12 5.743001 3.55032910
T12-1 12 21.974863 3.43095251
M67 12 28.863282 1.69070880
SSR345 12 44.236000 0.76826234
Tg50 12 55.627526 1.06820054
T12-3 12 73.950488 0.41372156
M68 12 90.821113 0.20482765
#重新排列每一條染色體的marker，找出最佳排列
> ripple(test01,chr=1)
204 total orders
obligXO
Initial 1 2 3 4 5 6 7 8 9 10 11 12 13 14 88
1 2 1 3 4 5 6 7 8 9 10 11 12 13 14 88
2 1 2 3 4 5 6 7 8 9 10 11 12 14 13 89
3 1 2 3 5 4 6 7 8 9 10 11 12 13 14 90
4 1 2 3 4 5 7 6 8 9 10 11 12 13 14 90
5 1 2 3 5 4 7 6 8 9 10 11 12 13 14 92
6 1 2 3 4 5 6 7 9 8 10 11 12 13 14 92
7 1 2 3 4 5 7 6 9 8 10 11 12 13 14 94
8 1 2 3 4 5 6 7 8 9 11 10 12 13 14 94
9 1 2 3 4 5 6 7 8 9 10 11 13 14 12 94
10 1 2 3 4 5 6 7 8 9 10 11 14 13 12 94
11 1 2 3 4 6 5 7 8 9 10 11 12 13 14 96
12 3 1 2 4 5 6 7 8 9 10 11 12 13 14 98
13 3 2 1 4 5 6 7 8 9 10 11 12 13 14 98
14 1 2 3 6 4 5 7 8 9 10 11 12 13 14 98
15 1 2 3 6 5 4 7 8 9 10 11 12 13 14 98
# MQM方法
> crossaug<-mqmaugment(test01)
> result<-mqmscan(crossaug)
> lodint(result,chr=2)
chr pos (cM) LOD pheno info LOD*info
T2-1 2 0 2.294957 0.9959819 2.285735
c2.loc10 2 10 3.173210 0.9376746 2.975438
c2.loc55 2 55 1.414865 0.9353935 1.323456
> lodint(result,chr=3)
chr pos (cM) LOD pheno info LOD*info
M11 3 0.00 0.02783887 1.0005814 0.02785506
c3.loc105 3 105.00 0.49768291 0.9880529 0.49173706
T3-2 3 129.06 0.10177260 1.0000000 0.10177260
#利用scanone的結果，去看QTL可能的區間
> sqt<-scanone(test01)
Warning message:
In scanone(test01) : First running calc.genoprob.
> lodint(sqt,chr=2)
chr pos lod
T2-1 2 0.00000 3.076447
OvtEx2 2 20.19844 5.360252
T2-2 2 46.71334 3.108854
#利用scanone的結果，利用貝氏法計算信來區間(Calculate an approximate Bayesian credible interval for a particular chromosome)
> bayesint(sqt,chr=2)
chr pos lod
M10 2 8.426426 4.227714
OvtEx2 2 20.198443 5.360252
SSR331 2 30.689005 4.650982
> bayesint(sqt,chr=12)
chr pos lod
M65 12 0.000000 2.996722
M66 12 5.743001 3.550329
T12-1 12 21.974863 3.430953
#Permutation test
#下載package>snow
>library(snow)
#permutation5000次，hk方法為regression of the phenotypes on the multipoint QTL genotype probabilities，verbose表示permutation的結果會
> a_perm<-scanone(test01,method="hk",n.perm=5000,verbose=T)
> ripple(test01,chr=1)
204 total orders
obligXO
Initial 1 2 3 4 5 6 7 8 9 10 11 12 13 14 88
1 2 1 3 4 5 6 7 8 9 10 11 12 13 14 88
2 1 2 3 4 5 6 7 8 9 10 11 12 14 13 89
3 1 2 3 5 4 6 7 8 9 10 11 12 13 14 90
4 1 2 3 4 5 7 6 8 9 10 11 12 13 14 90
5 1 2 3 5 4 7 6 8 9 10 11 12 13 14 92
6 1 2 3 4 5 6 7 9 8 10 11 12 13 14 92
7 1 2 3 4 5 7 6 9 8 10 11 12 13 14 94
8 1 2 3 4 5 6 7 8 9 11 10 12 13 14 94
9 1 2 3 4 5 6 7 8 9 10 11 13 14 12 94
10 1 2 3 4 5 6 7 8 9 10 11 14 13 12 94
11 1 2 3 4 6 5 7 8 9 10 11 12 13 14 96
12 3 1 2 4 5 6 7 8 9 10 11 12 13 14 98
13 3 2 1 4 5 6 7 8 9 10 11 12 13 14 98
14 1 2 3 6 4 5 7 8 9 10 11 12 13 14 98
15 1 2 3 6 5 4 7 8 9 10 11 12 13 14 98
# MQM方法
> crossaug<-mqmaugment(test01)
> result<-mqmscan(crossaug)
#顯示第二條染色體的QTL區間
> lodint(result,chr=2)
chr pos (cM) LOD pheno info LOD*info
T2-1 2 0 2.294957 0.9959819 2.285735
c2.loc10 2 10 3.173210 0.9376746 2.975438
c2.loc55 2 55 1.414865 0.9353935 1.323456
#Calculate a LOD support interval for a particular chromosome, using output from scanone.
> lodint(result,chr=3)
chr pos (cM) LOD pheno info LOD*info
M11 3 0.00 0.02783887 1.0005814 0.02785506
c3.loc105 3 105.00 0.49768291 0.9880529 0.49173706
T3-2 3 129.06 0.10177260 1.0000000 0.10177260
#利用scanone的結果，去看QTL可能的區間
> sqt<-scanone(test01)
Warning message:
In scanone(test01) : First running calc.genoprob.
> lodint(sqt,chr=2)
chr pos lod
T2-1 2 0.00000 3.076447
OvtEx2 2 20.19844 5.360252
T2-2 2 46.71334 3.108854
#利用scanone的結果，利用貝氏法計算信來區間(Calculate an approximate Bayesian credible interval for a particular chromosome)
> bayesint(sqt,chr=2)
M10 2 8.426426 4.227714
OvtEx2 2 20.198443 5.360252
SSR331 2 30.689005 4.650982
> bayesint(sqt,chr=12)
chr pos lod
M65 12 0.000000 2.996722
M66 12 5.743001 3.550329
T12-1 12 21.974863 3.430953
#Permutation test
#下載package>snow
>library(snow)
#permutation5000次，hk方法為regression of the phenotypes on the multipoint QTL genotype probabilities，verbose表示permutation的結果會
> a_threshold=summary(a_perm,alpha=c(0.1,0.05,0.01))
#calc.genoprob可以把marker之間每一個cM可能的互換率計算進去，故QTL計算的summary結果不是marker name而是染色體上的位置。
> sqt<-calc.genoprob(test01,step=1,map.function="kosambi")
> sqt<-calc.genoprob(test01,step=1,map.function="kosambi")
> out.fs<-scanone(sqt,method="hk")
> summary(out.fs,threshold=3)
chr pos lod
c2.loc17 2 17 5.71
c7.loc15 7 15 3.10
c12.loc14 12 14 4.03
> plot(out.fs)
> for(i in 1:3){abline(h=a_threshold[i],lty=3)}
> summary(out.fs,perm=a_perm,alpha=0.1)
chr pos lod
c2.loc17 2 17 5.71
c12.loc14 12 14 4.03
> summary(out.fs,perm=a_perm,alpha=0.05)
chr pos lod
c2.loc17 2 17 5.71
c12.loc14 12 14 4.03
> summary(out.fs,perm=a_perm,alpha=0.01)
chr pos lod
c2.loc17 2 17 5.71
#先把剛剛呈現最顯著LOD值的QTL加入
>B_qtl_1<-makeqtl(sqt,2,17,what="prob")
> plot(B_qtl_1)
#先把剛剛呈現最顯著LOD值的QTL加入
>B_qtl_1<-makeqtl(sqt,2,17,what="prob")
> plot(B_qtl_1)
> qc<-c(2,12)
> qp<-c(17,14)
> qtl<-makeqtl(sqt,qc,qp,what="prob")
> plot(qtl)
#看交感是否顯著
> addint(sqt,qtl=qtl,formula=y~Q1+Q2,method="hk")
Method: Haley-Knott regression
Model: normal phenotype
Model formula: y ~ Q1 + Q2
Add one pairwise interaction at a time table:
--------------------------------------------
df TypeIII SS LOD %var F value Pvalue(Chi2) Pvalue(F)
2@17.0:12@14.0 4 0.001834 0.2257 0.903 0.2114 0.904 0.931
Note:第二條17cM和第十二條14cM的交感不顯著
> out1<-fitqtl(sqt,qtl=qtl,method="hk")
> summary(out1)
fitqtl summary
Method: Haley-Knott regression
Model: normal phenotype
Number of observations : 46
Full model result
----------------------------------
Model formula: y ~ Q1 + Q2
df SS MS LOD %var Pvalue(Chi2) Pvalue(F)
Model 4 0.12103202 0.030258006 9.049756 59.58595 1.947393e-08 1.13497e-07
Error 41 0.08208972 0.002002188
Total 45 0.20312174
Drop one QTL at a time ANOVA table:
----------------------------------
df Type III SS LOD %var F value Pvalue(Chi2) Pvalue(F)
2@17.0 2 0.05364 5.023 26.41 13.396 0 3.33e-05 ***
12@14.0 2 0.03255 3.336 16.02 8.128 0 0.00106 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
沒有留言:
張貼留言