Single marker analysis

In this demonstration, I will show how to do single-marker analysis. This is test of association between a trait and a single marker.

Load data

load the phenotype and genotype data in the

load('barleyCAP.RData')

Select lines that were genotyped

pheno<- pheno[which(pheno$ID %in% row.names(geno)),]

Subset phenotype data

We will just work with one location and trait for this demonstration

pheno<- pheno[which(pheno$Site =='Genesee_2007'),]
pheno<- droplevels.data.frame(pheno[which(pheno$Trait == 'heading_date'),])

Subset the genotypic data

Want to get the marker data on the lines that have been phenotyped

geno2<- geno[pheno$ID,]

Look at the data

geno2[1:10, 1:10] #marker data

##             Mrk_3 Mrk_4 Mrk_5 Mrk_6 Mrk_7 Mrk_8 Mrk_9 Mrk_10 Mrk_11 Mrk_12
## 02MB-104       -1    -1     1    -1    -1     1     1      1     -1      1
## 03WA-109.7     -1     1     1     1     1    -1     1     -1     -1     -1
## 03WA-119.7      1    -1    -1    -1    -1    -1     1     -1      1     -1
## 03WA-124.1      1    -1     1    -1     1     1     1      1      1      1
## 03WA-129.18     1    -1     1    -1    -1    -1     1     -1     -1      1
## 03WA-158.12     1     1     1    -1    -1    -1     1      1     -1      1
## 04WA-101.34     1    -1     1    -1     1     1     1      1      1      1
## 04WA-101.45     1    -1     1    -1     1     1     1      1     -1      1
## 04WA-101.47     1    -1     1    -1     1     1     1      1      1      1
## 04WA-104.17    -1    -1     1    -1    -1    -1     1      1      1      1

head(pheno) #phenotype data

##      Exp     Env Year          ID        Trait         Site Value
## 8   CAP3 Genesee 2007    02MB-104 heading_date Genesee_2007    70
## 45  CAP1 Genesee 2007  03WA-109.7 heading_date Genesee_2007    60
## 84  CAP1 Genesee 2007  03WA-119.7 heading_date Genesee_2007    64
## 110 CAP1 Genesee 2007  03WA-124.1 heading_date Genesee_2007    63
## 123 CAP1 Genesee 2007 03WA-129.18 heading_date Genesee_2007    61
## 214 CAP1 Genesee 2007 03WA-158.12 heading_date Genesee_2007    62

Perform single marker analysis

Here we are determining the effect of marker i on the trait

i<-1
x<- geno2[,i] #
mod<- lm(pheno$Value~1+x) #regression model
eff<- mod$coefficients[2] #marker effect
smry<- summary(mod) #regression model summary
smry

## 
## Call:
## lm(formula = pheno$Value ~ 1 + x)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.704 -3.894  1.296  3.296  6.296 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  67.5837     0.3997 169.065   <2e-16 ***
## x            -0.1205     0.4040  -0.298    0.766    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.903 on 94 degrees of freedom
## Multiple R-squared:  0.0009462,  Adjusted R-squared:  -0.009682 
## F-statistic: 0.08903 on 1 and 94 DF,  p-value: 0.7661

Get the F-statsitic and, p-value

Here we are determining if marker i is significantlly associated with the trait

fstat<- smry$fstatistic[1] #f-statistic
df1<- smry$fstatistic[2] #numerator degrees of freedom
df2<- smry$fstatistic[3] #denominator degrees of freedom
pval<- pf(fstat, df1, df2, lower.tail = FALSE)
pval

##     value 
## 0.7660784