20190327_MakeCovariateFiles

workflowr

Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  analysis/20190327_MakeCovariateFiles.Rmd
    Untracked:  analysis/20190421_RegressOutRNASeqPCs.Rmd

Unstaged changes:
    Modified:   analysis/index.Rmd

After checking previous QC analyses of genotypes and RNA-seq data, this markdown will describe covariates used in gene expression modeling and the associated Rscript will output those covariates in the format required for MatrixEQTL. Specifically, the covariates that will be output are {NumGenotypePCs}, {NumRNASeqPCs}, sex, and perhaps RIN score.

The Rscript version of this R markdown is included in the snakemake workflow to generate the [{NumGenotypePCs}.{NumRNASeqPCs}.covariates] file using the knitr::purl() function.

library(plyr)
library(tidyverse)
library(knitr)
library(readxl)

Run this command to write this [.Rmd] into a [.R] that is included in the snakemake workflow for eQTL mapping:

# knitr::purl(input="../analysis/20190327_MakeCovariateFiles.Rmd", output="../analysis/20190327_MakeCovariateFiles.R")

# Use command line input to specify input and output if this is the Rscript version of this file (as opposed to Rmarkdown).
if(commandArgs()[4] == "--file=../../analysis/20190327_MakeCovariateFiles.R"){
  args <- commandArgs(trailingOnly = T)
  EmptyFamFilepath <- args[1]
  MetadataExcel <- args[2]
  CountFilepath <- args[3]
  NumRS.PCs <- args[4]
  GenotypePCs <- args[5]
  NumGenotype.PCs <- args[6]
  OutputFile <- args[7]
  print (args)

} else {
  CountFilepath <- '../output/CountTable.tpm.txt.gz'
  EmptyFamFilepath <- '../output/ForAssociationTesting.temp.fam'
  MetadataExcel <- '../data/Metadata.xlsx'
  GenotypePCs <- '../output/PopulationStructure/pca.eigenvec'
  NumGenotype.PCs <- 3
  NumRS.PCs <- 3
}

Read in data:

CountTable <- read.table(gzfile(CountFilepath), header=T, check.names=FALSE, row.names = 1)

# dimensions of count table
CountTable %>% dim()

[1] 50434    39

kable(head(CountTable))

EmptyFamFile <- read.table(EmptyFamFilepath, col.names=c("FID", "IID", "Father", "Mother", "SX", "Pheno"), stringsAsFactors = F) %>%
  select(-Pheno)

GenotypePCs.df <- read.table(GenotypePCs, header=T) %>% 
  select(-c("FID")) %>%
  rename_all(paste0, ".GT")
kable(head(GenotypePCs.df))

OtherMetadata <- as.data.frame(read_excel(MetadataExcel))
kable(head(OtherMetadata))

# RNA-seq PCA
pca_results <- CountTable %>%
  +0.1 %>%
  mutate(sumVar = rowSums(.)) %>%
  arrange(desc(sumVar)) %>%
  head(2500) %>%
  select(-sumVar) %>%
  log() %>%
  t() %>%
  prcomp(center=T, scale. = T)

RNASeqPCs.df <- pca_results$x %>% 
  as.data.frame %>% 
  rownames_to_column() %>%
  rename_all(paste0, ".RS")

Covariates <- EmptyFamFile %>%
  select(IID) %>%
  left_join(RNASeqPCs.df[0:NumRS.PCs+1], by=c("IID" = "rowname.RS")) %>%
  left_join(GenotypePCs.df[0:NumGenotype.PCs+1], by=c("IID" = "IID.GT")) %>%
  left_join(
    OtherMetadata %>% select(Individual.ID, SX),
    by = c("IID" = "Individual.ID")
  ) %>%
  mutate(SX = plyr::mapvalues(SX, from=c("M", "F"), to=c(0,1))) %>%
  t()

Warning: Column `IID`/`IID.GT` joining character vector and factor,
coercing into character vector

kable(head(Covariates))

Write table if this is the Rscript form of this Rmarkdown

if(commandArgs()[4] == "--file=../../analysis/20190327_MakeCovariateFiles.R"){
  write.table(Covariates, col.names = F, sep='\t', file=OutputFile, row.names=T, quote=F)
}

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
 [1] readxl_1.1.0    knitr_1.22      forcats_0.4.0   stringr_1.4.0  
 [5] dplyr_0.8.0.1   purrr_0.3.2     readr_1.3.1     tidyr_0.8.2    
 [9] tibble_2.1.1    ggplot2_3.1.0   tidyverse_1.2.1 plyr_1.8.4     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1       highr_0.8        cellranger_1.1.0 pillar_1.3.1    
 [5] compiler_3.5.1   git2r_0.24.0     workflowr_1.2.0  tools_3.5.1     
 [9] digest_0.6.18    lubridate_1.7.4  jsonlite_1.6     evaluate_0.13   
[13] nlme_3.1-137     gtable_0.3.0     lattice_0.20-38  pkgconfig_2.0.2 
[17] rlang_0.3.3      cli_1.1.0        rstudioapi_0.10  yaml_2.2.0      
[21] haven_2.1.0      xfun_0.6         withr_2.1.2      xml2_1.2.0      
[25] httr_1.4.0       hms_0.4.2        generics_0.0.2   fs_1.2.6        
[29] rprojroot_1.3-2  grid_3.5.1       tidyselect_0.2.5 glue_1.3.1      
[33] R6_2.4.0         rmarkdown_1.11   modelr_0.1.4     magrittr_1.5    
[37] backports_1.1.3  scales_1.0.0     htmltools_0.3.6  rvest_0.3.2     
[41] assertthat_0.2.1 colorspace_1.4-1 stringi_1.4.3    lazyeval_0.2.2  
[45] munsell_0.5.0    broom_0.5.1      crayon_1.3.4