20190327_MakeFamAndCovariateFiles

workflowr

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

Untracked files:
    Untracked:  analysis/20190321_Check-Kinship-And-PopulationStructure.Rmd
    Untracked:  analysis/20190325_MergingRNASeqLanes.Rmd
    Untracked:  analysis/20190326_Admixture.Rmd
    Untracked:  analysis/20190326_PCA.Rmd
    Untracked:  analysis/20190327_MakeFamAndCovariateFiles.Rmd
    Untracked:  analysis/20190327_MakeFamPhenotypeFile.Rmd
    Untracked:  docs/figure/20190321_Check-Kinship-And-PopulationStructure.Rmd/
    Untracked:  docs/figure/20190325_MergingRNASeqLanes.Rmd/
    Untracked:  docs/figure/20190326_Admixture.Rmd/
    Untracked:  docs/figure/20190326_PCA.Rmd/

Unstaged changes:
    Deleted:    ._workflowr.yml.swp
    Modified:   analysis/20190320_Check-RNAseq-PCs.Rmd
    Modified:   analysis/index.Rmd

this shows how I will output a plink [.fam] file containing phenotypes for eQTL association testing with gemma. The phenotypes (expresion matrix) which will be log-transformed TPM values, filtering for genes which have non-zero TPM values for all samples. The Rscript version of this R code is included in the snakemake workflow to generate the [.fam] file using the knitr::purl() function.

library(tidyverse)
library(knitr)

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_MakeFamPhenotypeFile.Rmd", output="../analysis/20190327_MakeFamPhenotypeFile.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_MakeFamPhenotypeFile.R"){
  args <- commandArgs(trailingOnly = T)
  CountFilepath <- args[1]
  EmptyFamFilepath <- args[2]
  PhenotypeOutFilepath <- args[3]
  PhenotypeListOutFilepath <- args[4]
  GenesBedFile <- args[5]
} else {
  CountFilepath <- '../output/CountTable.tpm.txt.gz'
  EmptyFamFilepath <- '../output/ForAssociationTesting.temp.fam'
  GenesBedFile <- '../data/cDNA.all.chromosomal.bed'
}

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)

# Will use GeneRegions to filter out non-autosomal genes
GeneChromosomes <- read.table(GenesBedFile, col.names=c("chromosome", "start", "stop", "gene", "score", "strand"), stringsAsFactors = F) %>%
  select(gene, chromosome)
kable(head(GeneChromosomes))

Here I will consider how to filter genes for association testing

GeneSet1 <- CountTable %>%
  rownames_to_column('gene') %>%
  mutate(Mean = rowMeans(select(., -gene))) %>%
  filter(Mean>1) %>%
  pull(gene)

GeneSet2 <- CountTable %>%
  rownames_to_column('gene') %>%
  filter_if(is.numeric, all_vars(.>0)) %>%
  pull(gene)

#Number genes left after filter method1
length(GeneSet1)

[1] 15696

#Number genes left after filter method2
length(GeneSet2)

[1] 18924

#Number genes left after intersection of both methods
length(intersect(GeneSet1, GeneSet2))

[1] 13852

Seems like a reasonable (and inclusive) way to filter genes (rather than eQTL testing all ~50000 transcripts, many of which contain mostly 0’s) is to just choose the genes with no 0 TPM values in any samples.

For now I will also filter out non-autosomal genes.

PhenotypesToOutput <- CountTable %>%
  rownames_to_column('gene') %>%
  filter_if(is.numeric, all_vars(.>0)) %>%
  merge(GeneChromosomes, by="gene", all.x=T) %>%
  filter(!chromosome %in% c("X", "Y", "MT")) %>%
  select(-chromosome) %>%
  column_to_rownames('gene') %>%
  log() %>%
  t()

row.names(PhenotypesToOutput) <- colnames(CountTable)


Output.df <- EmptyFamFile %>%
  merge(PhenotypesToOutput, all.x=T, by.x="IID", by.y=0) %>% as.tibble()

Warning: `as.tibble()` is deprecated, use `as_tibble()` (but mind the new semantics).
This warning is displayed once per session.

Output.df

# A tibble: 39 x 18,287
   IID   FID   Father Mother    SX ENSPTRT00000000… ENSPTRT00000000…
   <chr> <chr>  <int>  <int> <int>            <dbl>            <dbl>
 1 295   Pan_…      0      0     0           -0.503            0.709
 2 317   Pan_…      0      0     0           -0.933            0.598
 3 338   Pan_…      0      0     0           -1.08             1.37 
 4 389   Pan_…      0      0     0           -0.779            0.648
 5 438   Pan_…      0      0     0           -0.871            0.778
 6 456   Pan_…      0      0     0           -0.516            0.791
 7 462   Pan_…      0      0     0           -0.674            1.06 
 8 476   Pan_…      0      0     0           -1.09             0.668
 9 495   Pan_…      0      0     0            0.221            0.852
10 4x00… Pan_…      0      0     0           -0.797            1.36 
# … with 29 more rows, and 18,280 more variables:
#   ENSPTRT00000000050.4 <dbl>, ENSPTRT00000000057.4 <dbl>,
#   ENSPTRT00000000060.3 <dbl>, ENSPTRT00000000069.4 <dbl>,
#   ENSPTRT00000000071.5 <dbl>, ENSPTRT00000000083.6 <dbl>,
#   ENSPTRT00000000093.4 <dbl>, ENSPTRT00000000100.4 <dbl>,
#   ENSPTRT00000000105.6 <dbl>, ENSPTRT00000000106.3 <dbl>,
#   ENSPTRT00000000114.4 <dbl>, ENSPTRT00000000116.4 <dbl>,
#   ENSPTRT00000000118.6 <dbl>, ENSPTRT00000000144.6 <dbl>,
#   ENSPTRT00000000145.7 <dbl>, ENSPTRT00000000149.3 <dbl>,
#   ENSPTRT00000000152.5 <dbl>, ENSPTRT00000000172.4 <dbl>,
#   ENSPTRT00000000176.6 <dbl>, ENSPTRT00000000180.5 <dbl>,
#   ENSPTRT00000000188.6 <dbl>, ENSPTRT00000000207.3 <dbl>,
#   ENSPTRT00000000210.3 <dbl>, ENSPTRT00000000215.6 <dbl>,
#   ENSPTRT00000000216.4 <dbl>, ENSPTRT00000000219.6 <dbl>,
#   ENSPTRT00000000222.5 <dbl>, ENSPTRT00000000227.5 <dbl>,
#   ENSPTRT00000000229.3 <dbl>, ENSPTRT00000000234.5 <dbl>,
#   ENSPTRT00000000237.5 <dbl>, ENSPTRT00000000238.4 <dbl>,
#   ENSPTRT00000000243.6 <dbl>, ENSPTRT00000000247.4 <dbl>,
#   ENSPTRT00000000254.4 <dbl>, ENSPTRT00000000258.6 <dbl>,
#   ENSPTRT00000000260.4 <dbl>, ENSPTRT00000000270.4 <dbl>,
#   ENSPTRT00000000275.3 <dbl>, ENSPTRT00000000278.5 <dbl>,
#   ENSPTRT00000000279.6 <dbl>, ENSPTRT00000000298.3 <dbl>,
#   ENSPTRT00000000302.7 <dbl>, ENSPTRT00000000303.5 <dbl>,
#   ENSPTRT00000000306.5 <dbl>, ENSPTRT00000000310.5 <dbl>,
#   ENSPTRT00000000316.4 <dbl>, ENSPTRT00000000323.4 <dbl>,
#   ENSPTRT00000000324.4 <dbl>, ENSPTRT00000000329.4 <dbl>,
#   ENSPTRT00000000334.4 <dbl>, ENSPTRT00000000335.3 <dbl>,
#   ENSPTRT00000000336.5 <dbl>, ENSPTRT00000000339.4 <dbl>,
#   ENSPTRT00000000341.6 <dbl>, ENSPTRT00000000343.4 <dbl>,
#   ENSPTRT00000000348.3 <dbl>, ENSPTRT00000000349.4 <dbl>,
#   ENSPTRT00000000351.3 <dbl>, ENSPTRT00000000352.5 <dbl>,
#   ENSPTRT00000000354.6 <dbl>, ENSPTRT00000000355.4 <dbl>,
#   ENSPTRT00000000357.4 <dbl>, ENSPTRT00000000361.5 <dbl>,
#   ENSPTRT00000000369.4 <dbl>, ENSPTRT00000000395.4 <dbl>,
#   ENSPTRT00000000411.6 <dbl>, ENSPTRT00000000424.4 <dbl>,
#   ENSPTRT00000000431.5 <dbl>, ENSPTRT00000000435.3 <dbl>,
#   ENSPTRT00000000438.4 <dbl>, ENSPTRT00000000439.4 <dbl>,
#   ENSPTRT00000000441.5 <dbl>, ENSPTRT00000000444.5 <dbl>,
#   ENSPTRT00000000445.4 <dbl>, ENSPTRT00000000459.3 <dbl>,
#   ENSPTRT00000000462.4 <dbl>, ENSPTRT00000000464.4 <dbl>,
#   ENSPTRT00000000465.4 <dbl>, ENSPTRT00000000466.6 <dbl>,
#   ENSPTRT00000000494.5 <dbl>, ENSPTRT00000000499.4 <dbl>,
#   ENSPTRT00000000501.5 <dbl>, ENSPTRT00000000502.5 <dbl>,
#   ENSPTRT00000000510.3 <dbl>, ENSPTRT00000000524.5 <dbl>,
#   ENSPTRT00000000525.4 <dbl>, ENSPTRT00000000526.5 <dbl>,
#   ENSPTRT00000000532.4 <dbl>, ENSPTRT00000000535.3 <dbl>,
#   ENSPTRT00000000537.4 <dbl>, ENSPTRT00000000542.6 <dbl>,
#   ENSPTRT00000000548.4 <dbl>, ENSPTRT00000000555.5 <dbl>,
#   ENSPTRT00000000557.4 <dbl>, ENSPTRT00000000558.5 <dbl>,
#   ENSPTRT00000000568.4 <dbl>, ENSPTRT00000000569.4 <dbl>,
#   ENSPTRT00000000571.4 <dbl>, ENSPTRT00000000572.5 <dbl>, …

GeneList <- data.frame(GeneList=colnames(Output.df)[-1:-5])

if(commandArgs()[4] == "--file=../../analysis/20190327_MakeFamPhenotypeFile.R"){
  write.table(Output.df, col.names = F, sep='\t', file=PhenotypeOutFilepath, row.names=F, quote=F)
  write.table(GeneList, col.names = F, sep='\t', file=PhenotypeListOutFilepath, row.names=F, 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] knitr_1.22      forcats_0.4.0   stringr_1.4.0   dplyr_0.8.0.1  
 [5] purrr_0.3.2     readr_1.3.1     tidyr_0.8.2     tibble_2.1.1   
 [9] ggplot2_3.1.0   tidyverse_1.2.1

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