Last updated: 2020-03-30

Checks: 5 2

Knit directory: Blancetal-1/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.6.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: uncommitted changes

The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20200217)

The command set.seed(20200217) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute relative
~/Documents/Blancetal-1//analysis/ .

Repository version: 6b00f47

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 6b00f47. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 


Ignored files:
    Ignored:    analysis/.Rhistory

Untracked files:
    Untracked:  data/.~lock.gene_model_xref_v3.csv#
    Untracked:  data/quaint-results.rda
    Untracked:  data/siggenes.txt

Unstaged changes:
    Modified:   analysis/Identifying_quaint.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/Identifying_quaint.Rmd) and HTML (docs/Identifying_quaint.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 6b00f47 em 2020-03-27 stuff
html 6b00f47 em 2020-03-27 stuff
Rmd 640b45a em 2020-03-24 drought genes
html 640b45a em 2020-03-24 drought genes
Rmd 682f50e em 2020-03-23 stuff
html 682f50e em 2020-03-23 stuff
Rmd 450bded em 2020-03-19 more plots
html 450bded em 2020-03-19 more plots
Rmd 626d202 em 2020-03-19 num sig genes
Rmd 2688954 em 2020-03-18 merge conflict ugh
Rmd 8c52f3c Em 2020-03-17 stuff
Rmd dec95d3 Em 2020-03-04 pc figures
html dec95d3 Em 2020-03-04 pc figures
Rmd 993e4b7 Em 2020-03-03 mroe stuff
html 993e4b7 Em 2020-03-03 mroe stuff
Rmd a2e0aec Em 2020-02-27 drought gene analysis
html a2e0aec Em 2020-02-27 drought gene analysis
Rmd c018f1d Em 2020-02-26 added drought info
Rmd 227a4b6 Em 2020-02-26 added figure
html 227a4b6 Em 2020-02-26 added figure
Rmd a597c91 Em 2020-02-26 stuff
html a597c91 Em 2020-02-26 stuff

Intro

I’m rerunning the code to identify selected genes with the Quaint package

### function for testing all genes
testAllGenes <- function(myTissue){
# Read in mean-centered expression values
df1 <- read.table(paste("../data/Mean_centered_expression/",myTissue,".txt",sep=""))
geneNames = names(df1)
# Read in tissue specific kinship matrix 
myF <- read.table(paste('../data/Kinship_matrices/F_',myTissue,'.txt',sep=""))

## Get Eigen Values and Vectors 
myE <- eigen(myF)
E_vectors <- myE$vectors
E_values <- myE$values

## Testing for selection on first 5 PCs
myM = 1:5

## Using the last 1/2 of PCs to estimate Va
myL = 6:dim(myF)[1]

## test for selection on each gene
allGeneOutput <- sapply(1:ncol(df1), function(x){
myQpc = calcQpc(myZ = df1[,x], myU = E_vectors, myLambdas = E_values, myL = myL, myM = myM)
return(c(geneNames[x],myQpc))
})

return(allGeneOutput)
}

###run on all genes
alltissues = c('GRoot',"Kern","LMAD26","LMAN26", "L3Tip","GShoot","L3Base")
alltissueresults = lapply(alltissues, testAllGenes)
names(alltissueresults) <- alltissues



##look at sig results
sigresults = lapply(1:length(alltissues), function(i){
# extract the pvalues
pvals = matrix(unlist(alltissueresults[[i]][5,]), ncol=5, byrow=TRUE) #each row corresponds to a gene, columns are to PCs

##look at pvalue distributions
#par(mfrow=c(3,2), mar=c(4,4,1,1))
#sapply(1:5, function(x){
#hist(pvals[,x], border="white", col = "darkgray", main="", breaks = 20, xlab = paste('PC ',x,' ',alltissues[i],sep=""))
#})

## test fpr significance with the fdr
pfdr = data.frame(apply(pvals,2, function(x){p.adjust(x, method='fdr')}))

## how many are significant?
numsig <- apply(pfdr, 2, function(x){sum(x<0.1)})
numsig
})

###make a big image of how many sig results we have
sigTable = as.data.frame(matrix(unlist(sigresults), ncol=5, byrow=T))



names(sigTable) = c('PC1','PC2','PC3','PC4','PC5')
sigTable$nicename = c('Germinating root', 'Kernel','Adult Leaf Day', 'Adult Leaf Night', '3rd leaf tip', 'Germinating shoot','3rd leaf base')


sigLong = tidyr::gather(sigTable, 'variable','value', -nicename)
sigLong[sigLong == 0] <- NA

pl <- ggplot(data=sigLong,aes(x=variable,y=nicename)) + 
  geom_tile(aes(fill=value),color='black') + scale_fill_gradient(low = 'lightyellow', high = "#CC79A7", guide = FALSE, na.value = "white") + theme_bw() + xlab("\n") + ylab("Tissue") + 
  theme(axis.ticks=element_blank(),panel.border=element_blank(),panel.grid.major = element_blank(), axis.text.y = element_text(size=10,angle=0), axis.title.y = element_text(size=16),axis.title.x  = element_text(size=16),axis.text.x = element_text(angle = 0, hjust = 0.5,size=14)) + geom_text(aes(label=value),colour="grey15",size=5.5)

pl
Warning: Removed 20 rows containing missing values (geom_text).

###how many genes total?
sum(sigLong$value, na.rm=T) ##number of genes, but could have same gene for multiple tissues or PCs.
[1] 67
#unique genes
siggenes = unlist(lapply(1:length(alltissues), function(i){
# extract the pvalues
pvals = matrix(unlist(alltissueresults[[i]][5,]), ncol=5, byrow=TRUE) #each row corresponds to a gene, columns are to PCs
pfdr = data.frame(apply(pvals,2, function(x){p.adjust(x, method='fdr')})) ## test fpr significance with the fdr
pfdr$numsig <- apply(pfdr, 1, function(x){sum(x<0.1)})
pfdr$genes = unlist(alltissueresults[[i]][1,])
return(dplyr::filter(pfdr, numsig>0)$genes)
}))

length(unique(siggenes)) ##number of unique genes
[1] 60
###save output for other analyses
save(alltissueresults, sigTable, sigLong, file = "../data/quaint-results.rda")

Make some individual gene plots for leaf expression 26

##what gene to plot
pvals = matrix(unlist(alltissueresults[['LMAD26']][5,]), ncol=5, byrow=TRUE) #each row corresponds to a gene, columns are to 
pfdr = data.frame(apply(pvals,2, function(x){p.adjust(x, method='fdr')}))
pfdr$index = 1:nrow(pfdr)

##the sig LMAD26 one on PC 1
sigPC1 = dplyr::filter(pfdr, X1 < 0.1)




### function for making plots
makeGenePlot <- function(tissue, geneIndex, pc){

## read in population data
pop <- read.csv("../data/FlintGarciaTableS1_fixednew.csv", stringsAsFactors = F)

## read in expression data
exp <- read.table(paste("../data/Mean_centered_expression/",tissue,".txt", sep=""), stringsAsFactors = F)
expgene = data.frame(Inbred = row.names(exp), geneexp = exp[,geneIndex], stringsAsFactors = F)

## merge data together
pop_dat <- dplyr::inner_join(pop, expgene, by= "Inbred")

#get the eigenvalues
myF <- read.table(paste('../data/Kinship_matrices/F_',tissue,'.txt', sep=""))
myE = eigen(myF)
myPC = data.frame(pc = myE$vectors[,pc], stringsAsFactors = F)
pop_dat <- dplyr::bind_cols(pop_dat, myPC)

lambda <- myE$values[pc]


##calculate the CIs
generesults = alltissueresults[[tissue]][,geneIndex]
myVaest = var0(generesults$cml)
myCI = 1.96*sqrt(myVaest*lambda)

##gene name
geneName = names(exp)[geneIndex]

lR <- lm(pop_dat$geneexp ~ pop_dat$pc) 
coeff <- lR$coefficients[[2]]

if(tissue=="LMAD26"){mylab = c("mixed", "non-stiff stalk", "popcorn", "stiff stalk", "tropical")} else{mylab = c("mixed", "non-stiff stalk", "popcorn", "stiff stalk", "sweet", "tropical")} ##no sweets in LMAD26

col <- c('#E69F00', '#56B4E9', "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
pl1 <- ggplot(data=pop_dat, aes(x = pc, y= geneexp , color=Subpopulation)) + scale_colour_manual(values = col, labels=mylab) + xlab(paste("PC",pc)) +  ylab("Expression") + theme_bw() + theme(panel.grid.major = element_blank(), text = element_text(size=15), panel.grid.minor = element_blank(), axis.title.y = element_text(size=18), axis.title.x  = element_text(size=18), legend.position = "right", legend.title = element_text(size = 12), legend.text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + geom_point(size = 3.5) + geom_abline(slope = myCI, linetype = 2, col = '#56B4E9', size=1.2) + geom_abline(slope = coeff, size = 1.5)+  geom_abline(slope = -myCI, linetype = 2, col = '#56B4E9', size=1.2)

return(pl1)
}


myplot = makeGenePlot(tissue = 'LMAD26',geneIndex = sigPC1$index, pc = 1)

myplot

alltissueresults[['LMAD26']][1,sigPC1$index]##get gene name
[[1]]
[1] "GRMZM2G152686"
#make a plot for PC2
sigPC2 = dplyr::filter(pfdr, X2 < 0.1)
pc2plot = makeGenePlot(tissue="LMAD26",geneIndex = sigPC2$index[2], pc=2)
pc2plot

##make a plot for PC3
sigPC3 = dplyr::filter(pfdr, X3 < 0.1)
pc3plot = makeGenePlot(tissue="LMAD26",geneIndex = sigPC3$index[3], pc=3)
pc3plot

###make a plot for PC5
sigPC5 = dplyr::filter(pfdr, X5 < 0.1)
#sigPC5$gene = geneNames[sigPC5$index]
mgi = sigPC5$index[6]
alltissueresults[['LMAD26']][1,mgi]##get gene name
[[1]]
[1] "GRMZM2G069762"
pc5Plot = makeGenePlot(tissue='LMAD26',geneIndex=mgi, pc=5)

pc5Plot

##kernel tissue on PC 4
pvals = matrix(unlist(alltissueresults[['Kern']][4,]), ncol=5, byrow=TRUE) #each row corresponds to a gene, columns are to 
pfdr = data.frame(apply(pvals,2, function(x){p.adjust(x, method='fdr')}))
pfdr$index = 1:nrow(pfdr)
sigPC4 = dplyr::filter(pfdr, X4<0.1)

pc4Plot = makeGenePlot(tissue='Kern',geneIndex=sigPC4$index[1], pc=4)
pc4Plot

postscript("../figures/heatmap.eps",height=10,width=10,paper="special",horizontal=FALSE,colormodel="cymk", onefile=FALSE)

final <- ggarrange(pl,                                                 
          ggarrange(myplot, pc5Plot, ncol = 2, labels = c("B", "C"),common.legend = T, legend = "bottom"), 
          nrow = 2, 
          labels = "A"                                       
          )
Warning: Removed 20 rows containing missing values (geom_text).
final

dev.off()
png 
  2

Additional plots

mygene = "GRMZM2G093574"
#get gene index for kernel
kerngenes = unlist(alltissueresults[[2]][1,])
myindex = match(mygene, kerngenes)

#make a plot
explot = makeGenePlot('Kern',myindex, 2)
explot

##is this actually significant?
myps = matrix(unlist(alltissueresults[['Kern']][5,]), ncol=5, byrow=TRUE) #each row corresponds to 
myps[myindex,]
[1] 0.536121220 0.001042925 0.068229111 0.001874610 0.001091785
pfdr = data.frame(apply(myps,2, function(x){p.adjust(x, method='fdr')}))
pfdr[myindex,]
            X1        X2        X3        X4      X5
4039 0.9621721 0.5319451 0.9997441 0.9999446 0.99998

Look for interesting genes

allGenes = lapply(alltissues, function(mytissue){
myresults = alltissueresults[[match(mytissue, alltissues)]]
mygenes = unlist(myresults[1,])
myps = do.call(rbind,myresults[5,])
outdf = data.frame(mygenes, myps, stringsAsFactors = F)
names(outdf) = c('gene','PC1','PC2','PC3','PC4','PC5')
outdflong = tidyr::gather(outdf, 'PC','pval', PC1:PC5)
myfdrs = apply(myps, 2, p.adjust, method='fdr') ##get the fdrs
outdf2 = data.frame(mygenes, myfdrs, stringsAsFactors = F)
names(outdf2) = c('gene','PC1','PC2','PC3','PC4','PC5')
outdflong2 = tidyr::gather(outdf2, 'PC','fdr', PC1:PC5)
finaldf = dplyr::left_join(outdflong, outdflong2)
finaldf$tissue = mytissue
return(finaldf)
})
Joining, by = c("gene", "PC")
Joining, by = c("gene", "PC")
Joining, by = c("gene", "PC")
Joining, by = c("gene", "PC")
Joining, by = c("gene", "PC")
Joining, by = c("gene", "PC")
Joining, by = c("gene", "PC")
allgenes = do.call(rbind, allGenes)

siggenes = dplyr::filter(allgenes, fdr < 0.1)

##any with more than one association?
siggenes[duplicated(siggenes$gene),]
            gene  PC         pval          fdr tissue
15 GRMZM2G152686 PC2 8.588522e-06 0.0525904162 LMAD26
47 GRMZM2G152686 PC1 1.801337e-08 0.0001519428 LMAN26
53 GRMZM5G829894 PC3 1.882318e-05 0.0645920908 LMAN26
59 GRMZM2G391936 PC5 3.031756e-06 0.0127864312 LMAN26
61 GRMZM5G800586 PC5 6.661297e-05 0.0936467291 LMAN26
64 GRMZM2G053801 PC2 1.469630e-06 0.0149828809 GShoot
66 GRMZM2G093574 PC2 2.717956e-06 0.0314059793 L3Base
write.table(allgenes, file = '../data/siggenes.txt', quote=F, row.names=F, col.names=F)

https://www.maizegdb.org/gene_center/gene/GRMZM2G152686 https://www.maizegdb.org/gene_center/gene/GRMZM2G391936

Deought tolerance!! https://www.maizegdb.org/gene_center/gene/GRMZM5G800586


sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.6 LTS

Matrix products: default
BLAS:   /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] tidyr_1.0.2       quaint_0.0.0.9000 ggpubr_0.2.5      magrittr_1.5     
[5] reshape2_1.4.3    ggplot2_3.3.0    

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4       plyr_1.8.6       compiler_3.6.3   pillar_1.4.3    
 [5] later_1.0.0      git2r_0.26.1     highr_0.8        workflowr_1.6.1 
 [9] tools_3.6.3      digest_0.6.25    evaluate_0.14    lifecycle_0.2.0 
[13] tibble_2.1.3     gtable_0.3.0     pkgconfig_2.0.3  rlang_0.4.5     
[17] rstudioapi_0.11  yaml_2.2.1       xfun_0.12        gridExtra_2.3   
[21] withr_2.1.2      stringr_1.4.0    dplyr_0.8.5      knitr_1.28      
[25] vctrs_0.2.4      fs_1.3.2         cowplot_1.0.0    rprojroot_1.3-2 
[29] grid_3.6.3       tidyselect_1.0.0 glue_1.3.2       R6_2.4.1        
[33] rmarkdown_2.1    farver_2.0.3     purrr_0.3.3      whisker_0.4     
[37] ellipsis_0.3.0   backports_1.1.5  scales_1.1.0     promises_1.1.0  
[41] htmltools_0.4.0  assertthat_0.2.1 colorspace_1.4-1 ggsignif_0.6.0  
[45] httpuv_1.5.2     labeling_0.3     stringi_1.4.6    munsell_0.5.0   
[49] crayon_1.3.4