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.3.1   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.1.1

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4  haven_1.1.2       lattice_0.20-35  
 [4] colorspace_1.3-2  htmltools_0.3.6   yaml_2.2.0       
 [7] rlang_0.2.2       R.oo_1.22.0       pillar_1.3.0     
[10] glue_1.3.0        withr_2.1.2       R.utils_2.7.0    
[13] modelr_0.1.2      readxl_1.1.0      bindr_0.1.1      
[16] plyr_1.8.4        munsell_0.5.0     gtable_0.2.0     
[19] cellranger_1.1.0  rvest_0.3.2       R.methodsS3_1.7.1
[22] evaluate_0.11     knitr_1.20        broom_0.5.0      
[25] Rcpp_0.12.19      scales_1.0.0      backports_1.1.2  
[28] jsonlite_1.5      hms_0.4.2         digest_0.6.17    
[31] stringi_1.2.4     grid_3.5.1        rprojroot_1.3-2  
[34] cli_1.0.1         tools_3.5.1       magrittr_1.5     
[37] lazyeval_0.2.1    crayon_1.3.4      whisker_0.3-2    
[40] pkgconfig_2.0.2   xml2_1.2.0        lubridate_1.7.4  
[43] assertthat_0.2.0  rmarkdown_1.10    httr_1.3.1       
[46] rstudioapi_0.8    R6_2.3.0          nlme_3.1-137     
[49] git2r_0.23.0      compiler_3.5.1