R/qtl 找QTL與Confidence interval方法分享

(一)資料整理

橫列依序為:1:marker name、2:位於第幾條染色體、3:先用MapDisto或Mapmaker計算遺傳距離、4~50:每一個marker的genotype data
直行依序為: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的找尋。

> 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.07644668
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

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)

#顯示第二條染色體的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) 

              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)

> 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")

> 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)

> 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 

R: Dunnett test 方法

假設共有4組試驗, 其中F1(LA1312xLA1543)為對照組, 另有3 組實驗組(親本LA1312,LA1543,LA1589)。每組試驗各有5重複。

1.安裝package-> multcomp

#20個番茄果實長寬比的觀測值
>y=c(1.233754545,1.1738625,1.21531,1.3708,1.251022222,1.464914286,1.4413,1.266463636,1.492966667,1.30975,1.302936364,1.296875,1.2109,1.286966667,1.2878,1.013746667,0.986116667,1.008646154,1.0426,1.031922222)

#將其分為四組，共五重複
> trt=gl(4,5,label=c("F1","LA1312","LA1543","LA1589"))

#打開package，沒有打開就做不了，所以一定要用library打開。
> library(multcomp)
Loading required package: mvtnorm
Loading required package: survival
Loading required package: splines

#執行ANOVA，發現其有顯著差異，再利用Dunnett 去看是哪幾個跟對照組有差異。
> out=aov(y~trt)
> summary(out)
            Df Sum Sq Mean Sq F value   Pr(>F)    
trt          3 0.3764 0.12548   28.75 1.11e-06 ***
Residuals   16 0.0698 0.00436                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

#glht 來做單一對比與直交對比
> glht(out,linfct=mcp(trt="Dunnett"))

         General Linear Hypotheses


Multiple Comparisons of Means: Dunnett Contrasts




Linear Hypotheses:
                 Estimate
LA1312 - F1 == 0  0.14613
LA1543 - F1 == 0  0.02815
LA1589 - F1 == 0 -0.23234

> dunnett.out=glht(out,linfct=mcp(trt="Dunnett"))
> summary(dunnett.out)

         Simultaneous Tests for General Linear Hypotheses


Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = y ~ trt)


Linear Hypotheses:
                 Estimate Std. Error t value Pr(>|t|)    
LA1312 - F1 == 0  0.14613    0.04178   3.497  0.00837 ** 
LA1543 - F1 == 0  0.02815    0.04178   0.674  0.84038    
LA1589 - F1 == 0 -0.23234    0.04178  -5.561  < 0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- single-step method)


結論:可知F1(LA1543xLA1312)與LA1543形狀相似，但與LA1312和LA1589形狀有顯著不同。

土壤學

磷肥
土壤鹼性性:磷酸鈣
土壤酸性:磷酸鐵、磷酸鋁(磷被固定)
移動性低，所以要施用在根系旁邊，條撥予點施。
與氮素施用與有機質施用，土壤過鹼:硫酸鈣(石膏)，土壤過酸:添加石灰

氮肥:
尿素中氮的含量:46%
硫酸氨氮的含量:21%

重量含水量: 風乾土-烘乾土/烘乾土
總體積密度:風乾土重/風乾土體積
孔隙率:1-(土壤密度/粒子密度)


微生物分解有機質流程:(查!)
H2S>H2SO4
NH3>NH4

基肥>追肥>禮肥(多年生植物，收穫後以補足其果實消耗之養分)

砂土
缺點:保水保肥能力差、有機質少
改進:加入保水保肥的土壤

黏土:
改進:深耕把土壤打碎或加入砂土

土壤膠體性質:
無機膠體:
金屬、結晶矽酸鹽、非結晶矽酸鹽
有機膠體:
有機質

五種氮肥:
硫酸銨:生理酸性，不能連續使用以免土壤酸化，可加石灰以中和其酸性(不可施用於石灰質土壤中，以免氮素流失)
尿素:易流失，易溶於水
銨態氮:緩效肥
氫氮化鈣:對新芽有毒害，為鹼性

P-R=I=E+T+D+W
w:土中水分改變的量，通常為0
D:土壤剖面滲流出的量
E:土表蒸發
T:葉面蒸散
P:土壤水分總體積(天然降雨或是灌溉水)
R:水分逕流量

有機質過多，會使土壤酸化，
可添加有益微生物於土壤中，分解有機物質的細菌，抗生素細菌，輪作栽培，調整土壤pH值(6-7)，適量施用農藥，保持微生物生存空。

土壤好壞要如何判定:

1.土壤有芳香味，黑棕色
2.種子發芽法判定，
發芽試驗=堆肥抽出液中發芽率/去離子水中發芽率。
3.塑膠袋法，未成熟，氣體會膨大，若以為完全腐熟，則較少
4.蚯蚓法，如果蚯蚓往下鑽，則為腐熟堆肥
5.測pH值>7為好的堆肥，速測土壤pH值，土壤體積:水= 1:5
6.測全氮、全鉀含量

黃棕色:鐵礦
灰色:若欲吐，土壤進水
紅色:鐵的氧化物，赤鐵礦，偏紅色土壤化育層較高，受地下水位影響，鐵氧化物在水位變動地方累積。
黑色:有機質和猛氧化物的結合，地下水變動地方
白色:河床或石灰質土壤，海邊鹽類結晶，
青灰色:頁岩、泥岩顏色，排水不良(藍色)->還原態

台灣農田有機質含量:1%~3% >3%

複合肥:多用於基肥，因其個元素含量固定較難配和各生長期。
單質肥:多用於追肥 。

研究策略

(一)要如何證明mRNA有binding到DNA上，或microRNA對mRNA有影響?

1.Wobble mutant: 改變mRNA構型，但不改變protein，利用third codon的mutant。

例子: 看micro RNA 是否對RNA有影響?
方法: 改變第三個密碼子，使蛋白質相同，但mricoRNA不能切。

(二)Southern的應用:

1.檢視某一gene的family大小，是單組或多組，如果是多組就要考慮其他基因family的影響，是否有duplication?或有幾個copy?

Epigenetic 後生遺傳學

Epigenetic:

定義:
相同DNA背景，卻有不同的mRNA表現，可能是histion(在Histion tail的修飾作用)
,micro RNA, siRNA,或DNA甲基化的原因(使染色質緊緊纏繞，導致基因不表現)。

維基百科說明，此為研究DNA序列中未包含的基因調控訊息如何傳遞到下一代。

Genome-Wide Association Mapping

Association mapping also known as "linkage disequilibrium mapping", is a method of mapping quantitative trait loci (QTLs) that takes advantage of historic linkage disequilibrium to link phenotypes(observable characteristics) to genotypes (the genetic constitution of organisms).

Genomic selection in Plant Breeding : 
(出處 Advance in Agronomy Volume 110, 2011, Pages 77–123)




“Genomic selection,” the ability to select for even complex, quantitative traits based on marker data alone, has arisen from the conjunction of new high-throughput marker technologies and new statistical methods needed to analyze the data. This review surveys what is known about these technologies, with sections on population and quantitative genetic background, DNA marker development, statistical methods, reported accuracies of genomic selection (GS) predictions, prediction of nonadditive genetic effects, prediction in the presence of subpopulation structure, and impacts of GS on long-term gain. GS works by estimating the effects of many loci spread across the genome. Marker and observation numbers therefore need to scale with the genetic map length in Morgans and with the effective population size of the population under GS. For typical crops, the requirements range from at least 200 to at most 10,000 markers and observations. With that baseline, GS can greatly accelerate the breeding cycle while also using marker information to maintain genetic diversity and potentially prolong gain beyond what is possible with phenotypic selection. With the costs of marker technologies continuing to decline and the statistical methods becoming more routine, the results reviewed here suggest that GS will play a large role in the plant breeding of the future. Our summary and interpretation should prove useful to breeders as they assess the value of GS in the context of their populations and resources.