Exploritory analysis prior to running ash

Apa and Protein

Upload data:

I want the filtered peak counts. I need to filter the counts file for the total fraction based on the filtered peaks.

Total counts: /project2/gilad/briana/threeprimeseq/data/filtPeakOppstrand_cov_noMP/filtered_APApeaks_merged_allchrom_refseqGenes.TranscriptNoMP_sm_quant.Total_fixed.fc

okPeaks: /project2/gilad/briana/threeprimeseq/data/PeakUsage_noMP/filtered_APApeaks_merged_allchrom_refseqGenes.TranscriptNoMP_sm_quant.Total_fixed.pheno.5percPeaks.txt

filterTotalCounts_noMP_5percCov.py

totalokPeaks5perc_file="/project2/gilad/briana/threeprimeseq/data/PeakUsage_noMP/filtered_APApeaks_merged_allchrom_refseqGenes.TranscriptNoMP_sm_quant.Total_fixed.pheno.5percPeaks.txt"
countFile=open("/project2/gilad/briana/threeprimeseq/data/filtPeakOppstrand_cov_noMP/filtered_APApeaks_merged_allchrom_refseqGenes.TranscriptNoMP_sm_quant.Total_fixed.fc","r")

outFile=open("/project2/gilad/briana/threeprimeseq/data/filtPeakOppstrand_cov_noMP/filtered_APApeaks_merged_allchrom_refseqGenes.TranscriptNoMP_sm_quant.Total_5percCov_fixed.fc", "w")

allPeakOk={}
for ln in open(totalokPeaks5perc_file,"r"):
    peakname=ln.split()[5]
    peaknum=peakname[4:]
    if peaknum not in allPeakOk.keys():
        allPeakOk[peaknum]=""
    
for i,ln in enumerate(countFile):
    if i==1:
      outFile.write(ln)
    if i>1:
      ID=ln.split()[0]
      peak=ID.split(":")[0]
      peaknum=peak[4:]
      if peaknum in allPeakOk.keys():
          outFile.write(ln)
outFile.close()

total_Cov=read.table("../data/PeakCounts_noMP_5perc/filtered_APApeaks_merged_allchrom_refseqGenes.TranscriptNoMP_sm_quant.Total_5percCov_fixed.fc", stringsAsFactors = F,header = T) %>%  separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "Name"), sep=":")
total_genes=total_Cov %>% select(Name) %>% arrange(Name) %>% unique()

Gene names:

geneNames=read.table("../data/ensemble_to_genename.txt",sep="\t", header=T, stringsAsFactors = F, col.names=c("ID", "Name", "Source"))

prot=read.table("../data/mol_pheno/fastqtl_qqnorm_prot.fixed.noChr.txt",header=T,stringsAsFactors = F) %>% inner_join(geneNames, by="ID")

Keep the protein genes in APA:

prot_inAPA=prot %>% semi_join(total_genes, by="Name")

This shows we have 4209 genes with data for both. Now I can back filter the total peaks for the genes in prot_inAPA

total_Cov_wProt= total_Cov %>% semi_join(prot_inAPA,by="Name")

Need to give Stephens lab: unadjusted R-squared, and n and p for each protein, where p is the number of apa s that you are using in the regression and n is the number of samples?

To do this I need to get the overlapping individuals:

protInd=colnames(prot)[5:(dim(prot)[2]-2)]

ApaInd=c()
for (i in colnames(total_Cov)[12:ncol(total_Cov)]){
  num=substr(i,2,6)
  name=paste("NA", num, sep="")
  ApaInd=c(ApaInd, name)
}


IndBoth=intersect(protInd,ApaInd)
length(IndBoth)

[1] 29

I have 29 individuals in common for these.

I need to make a matrix for each gene. It will have a row for each commmon individual. A column for the protein, and a column for each assocaited peaks. After I have this I will be able to get the R2 value.

First create a function.

get_R2=function(gene, Cov=total_Cov, prot=prot_inAPA, apaName=ApaInd){
  gene_un= enexpr(gene)
  #deal with APA
  genePeaks=total_Cov %>% filter(Name==!!gene_un)
  n=nrow(genePeaks)
  drop_col=c('chr','Chr', 'start','end','strand','Name', 'Start','End','Strand','Length')
  genePeaks_sm= genePeaks %>% select(-one_of(drop_col))
  colnames(genePeaks_sm)=c("peak", ApaInd) 
  genePeakM=genePeaks_sm %>% column_to_rownames(var="peak") %>% t()
  genePeakDF=as.data.frame(genePeakM) %>% rownames_to_column(var="Ind")
  #deal with prot
  drop_col_prot= c("Chr", "start", "end", "ID", "Name", "Source")
  geneProt=prot %>% filter(Name==!!gene_un) %>% select(-one_of(drop_col_prot)) %>% t()
  colnames(geneProt)="prot"
  #print(dim(geneProt))
  geneProt_df=as.data.frame(geneProt) %>% rownames_to_column(var="Ind") %>% drop_na(prot)
  #print(geneProt_df)
  both=geneProt_df %>% inner_join(genePeakDF,by="Ind")
  num=seq(1,n)
  base="summary(lm(both$prot~"
  for (i in 3:dim(both)[2]){
    base=paste(base, "+both[,",i,"]",sep="")
  }
  code=paste(base, "))$r.squared", sep="")
  r2=eval(parse(text=code))
  final=c(gene, r2,nrow(both),n)
  return(final)
}

Run this on all genes:

final_matrix=matrix(c("gene","r2","n","p"),1,4)
for (i in prot_inAPA$Name){
 final_matrix= rbind(final_matrix,get_R2(i))
}

Make this a dataframe:

final_df=as.data.frame(final_matrix)

colnames(final_df)=as.character(unlist(final_df[1,]))

final_df <- final_df[-1 ,]

save(final_df,file="../data/protAndAPAlmRes.Rda")

When the stephens lab ran ASH on this, all R2 shrunk to zero. ##Protein and expression

I want to look at the the protein ~ expression model. This is easier because there is always only 1 expression level per gene.

rna=read.table("../data/mol_pheno/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt",header=T,stringsAsFactors = F) %>%  separate(ID, into=c("ID", "ver"), sep ="[.]") %>%  inner_join(geneNames, by="ID")

I want to filter this by genes we have proteinn for.

rnaandProt= rna %>% semi_join(prot, by="Name")

These are the gene I want to run the analsis on. I will make a similar function to run the linear model.

get_R2_Protexp=function(gene, exp=rnaandProt, protein=prot){
  #gene_un= enexpr(gene)
  exp_gene=exp %>% filter(Name ==gene)
  #print(exp_gene)
  drop_col_exp= c("Chr", "start", "end", "ID", "ver", "Name", "Source")
  exp_gene_sm= exp_gene %>% select(-one_of(drop_col_exp)) %>% t()
  colnames(exp_gene_sm)="Expression"
  exp_gene_df=as.data.frame(exp_gene_sm) %>% rownames_to_column(var="Ind")
  #print(exp_gene_df)
  drop_col_p= c("Chr", "start", "end", "ID", "Name", "Source")
  prot_gene= protein %>% filter(Name ==gene)%>% select(-one_of(drop_col_p)) %>% t()
  colnames(prot_gene)= "Protein"
  prot_gene_df= as.data.frame(prot_gene) %>% rownames_to_column(var="Ind") %>% drop_na(Protein)
  #print(prot_gene_df)
  both=prot_gene_df %>% inner_join(exp_gene_df, by="Ind")
  #print(both)
  r2=summary(lm(both$Protein ~both$Expression))$r.squared
  #print(r2)
  final=c(gene, r2, nrow(both))
}

test=get_R2_Protexp(gene="ISG15")

Run on all genes in rnaandProt

final_matrix_protExp=matrix(c("gene","r2","n"),1,3)
for (i in rnaandProt$Name){
 final_matrix_protExp= rbind(final_matrix_protExp,get_R2_Protexp(i))
}

Fix as df and save

final_df_protExp=as.data.frame(final_matrix_protExp)

colnames(final_df_protExp)=as.character(unlist(final_df_protExp[1,]))

final_df_protExp <- final_df_protExp[-1 ,]

save(final_df_protExp,file="../data/protAndExpressionlmRes.Rda")

Protein, APA, Exp

I need to subset the protein for genes in apa and expr.

prot_inAPAandExp=prot %>% semi_join(total_genes, by="Name") %>% semi_join(rna,by="Name")

get_R2_full=function(gene, Cov=total_Cov, prot=prot_inAPAandExp, apaName=ApaInd,exp=rna){
  gene_un= enexpr(gene)
  #deal with APA
  genePeaks=total_Cov %>% filter(Name==!!gene_un)
  n=nrow(genePeaks)
  drop_col=c('chr','Chr', 'start','end','strand','Name', 'Start','End','Strand','Length')
  genePeaks_sm= genePeaks %>% select(-one_of(drop_col))
  colnames(genePeaks_sm)=c("peak", ApaInd) 
  genePeakM=genePeaks_sm %>% column_to_rownames(var="peak") %>% t()
  genePeakDF=as.data.frame(genePeakM) %>% rownames_to_column(var="Ind")
  #deal with prot
  drop_col_prot= c("Chr", "start", "end", "ID", "Name", "Source")
  geneProt=prot %>% filter(Name==!!gene_un) %>% select(-one_of(drop_col_prot)) %>% t()
  colnames(geneProt)="prot"
  #print(dim(geneProt))
  geneProt_df=as.data.frame(geneProt) %>% rownames_to_column(var="Ind") %>% drop_na(prot)
  #print(geneProt_df)
  #deal with expr
  exp_gene=exp %>% filter(Name ==gene)
  drop_col_exp= c("Chr", "start", "end", "ID", "ver", "Name", "Source")
  exp_gene_sm= exp_gene %>% select(-one_of(drop_col_exp)) %>% t()
  colnames(exp_gene_sm)="Expression"
  exp_gene_df=as.data.frame(exp_gene_sm) %>% rownames_to_column(var="Ind")
  #make full model  
  both=geneProt_df %>% inner_join(exp_gene_df, by="Ind") %>% inner_join(genePeakDF,by="Ind")
  num=seq(1,n)
  base="summary(lm(both$prot~ both$Expression"
  for (i in 4:dim(both)[2]){
    base=paste(base, "+both[,",i,"]",sep="")
  }
  code=paste(base, "))$r.squared", sep="")
  #print(code)
  r2=eval(parse(text=code))
  final=c(gene, r2,nrow(both),n)
  return(final)
}

test_full=get_R2_full("ISG15")

Run this on all genes:

final_matrix_full=matrix(c("gene","r2","n","p"),1,4)
for (i in prot_inAPAandExp$Name){
 final_matrix_full= rbind(final_matrix_full,get_R2_full(i))
}

Make this a dataframe:

final_df_full=as.data.frame(final_matrix_full)

colnames(final_df_full)=as.character(unlist(final_df_full[1,]))

final_df_full <- final_df_full[-1 ,]

save(final_df_full,file="../data/protAndAPAAndExplmRes.Rda")

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

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] bindrcpp_0.2.2  cowplot_0.9.3   reshape2_1.4.3  forcats_0.3.0  
 [5] stringr_1.4.0   dplyr_0.7.6     purrr_0.2.5     readr_1.1.1    
 [9] tidyr_0.8.1     tibble_1.4.2    ggplot2_3.0.0   tidyverse_1.2.1
[13] workflowr_1.2.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.19     cellranger_1.1.0 plyr_1.8.4       compiler_3.5.1  
 [5] pillar_1.3.0     git2r_0.24.0     bindr_0.1.1      tools_3.5.1     
 [9] digest_0.6.17    lubridate_1.7.4  jsonlite_1.6     evaluate_0.13   
[13] nlme_3.1-137     gtable_0.2.0     lattice_0.20-35  pkgconfig_2.0.2 
[17] rlang_0.2.2      cli_1.0.1        rstudioapi_0.9.0 yaml_2.2.0      
[21] haven_1.1.2      withr_2.1.2      xml2_1.2.0       httr_1.3.1      
[25] knitr_1.20       hms_0.4.2        fs_1.2.6         rprojroot_1.3-2 
[29] grid_3.5.1       tidyselect_0.2.4 glue_1.3.0       R6_2.3.0        
[33] readxl_1.1.0     rmarkdown_1.11   modelr_0.1.2     magrittr_1.5    
[37] whisker_0.3-2    backports_1.1.2  scales_1.0.0     htmltools_0.3.6 
[41] rvest_0.3.2      assertthat_0.2.0 colorspace_1.3-2 stringi_1.2.4   
[45] lazyeval_0.2.1   munsell_0.5.0    broom_0.5.0      crayon_1.3.4