Understand Peaks

Variation in peaks

First I will run PCA on the peak coverage. I will run this seperatly for the total and nuclear fractions. I do not expect large amount of separation.

I will use the peak coverage data before the ratios are created for leafcutter. These files were created using feature counts on the filtered peaks. At this point the peaks have been mapped to the closest refseq transcript on the opposite strand.

Relevant file:
* /project2/gilad/briana/threeprimeseq/data/filtPeakOppstrand_cov/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc

• /project2/gilad/briana/threeprimeseq/data/filtPeakOppstrand_cov/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc

These files are in /Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/PeakCounts on my computer.

library(tidyverse)

── Attaching packages ───────────────────────────────────────────────────────── tidyverse 1.2.1 ──

✔ ggplot2 3.0.0     ✔ purrr   0.2.5
✔ tibble  1.4.2     ✔ dplyr   0.7.6
✔ tidyr   0.8.1     ✔ stringr 1.3.1
✔ readr   1.1.1     ✔ forcats 0.3.0

── Conflicts ──────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(workflowr)

This is workflowr version 1.1.1
Run ?workflowr for help getting started

library(cowplot)

Attaching package: 'cowplot'

The following object is masked from 'package:ggplot2':

    ggsave

library(reshape2)

Attaching package: 'reshape2'

The following object is masked from 'package:tidyr':

    smiths

library(devtools)
library(tximport)

Load data:

#only keep the counts 
total_Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,7:45]
nuclear_Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,7:45]

ggplot(total_Cov, aes(x=log10(X18486_T))) + geom_density()

Warning: Removed 233009 rows containing non-finite values (stat_density).

pca_tot_peak=prcomp(total_Cov, center=T,scale=T)
summary(pca_tot_peak)

Importance of components:
                          PC1     PC2     PC3     PC4     PC5    PC6
Standard deviation     5.9010 1.30000 0.81376 0.75658 0.47993 0.4501
Proportion of Variance 0.8929 0.04333 0.01698 0.01468 0.00591 0.0052
Cumulative Proportion  0.8929 0.93621 0.95319 0.96787 0.97378 0.9790
                           PC7     PC8     PC9    PC10    PC11    PC12
Standard deviation     0.42896 0.32313 0.30419 0.27984 0.23427 0.19916
Proportion of Variance 0.00472 0.00268 0.00237 0.00201 0.00141 0.00102
Cumulative Proportion  0.98369 0.98637 0.98874 0.99075 0.99216 0.99317
                          PC13    PC14    PC15    PC16    PC17   PC18
Standard deviation     0.18883 0.15913 0.15127 0.14309 0.12758 0.1254
Proportion of Variance 0.00091 0.00065 0.00059 0.00053 0.00042 0.0004
Cumulative Proportion  0.99409 0.99474 0.99532 0.99585 0.99626 0.9967
                          PC19    PC20    PC21    PC22    PC23    PC24
Standard deviation     0.12328 0.11035 0.10707 0.09979 0.09530 0.08797
Proportion of Variance 0.00039 0.00031 0.00029 0.00026 0.00023 0.00020
Cumulative Proportion  0.99706 0.99737 0.99766 0.99792 0.99815 0.99835
                          PC25    PC26    PC27    PC28    PC29    PC30
Standard deviation     0.08576 0.08086 0.07902 0.07535 0.07454 0.06907
Proportion of Variance 0.00019 0.00017 0.00016 0.00015 0.00014 0.00012
Cumulative Proportion  0.99854 0.99871 0.99887 0.99901 0.99916 0.99928
                          PC31    PC32    PC33    PC34    PC35    PC36
Standard deviation     0.06717 0.06441 0.06201 0.05666 0.05415 0.05261
Proportion of Variance 0.00012 0.00011 0.00010 0.00008 0.00008 0.00007
Cumulative Proportion  0.99939 0.99950 0.99960 0.99968 0.99976 0.99983
                          PC37    PC38    PC39
Standard deviation     0.05128 0.04839 0.04237
Proportion of Variance 0.00007 0.00006 0.00005
Cumulative Proportion  0.99989 0.99995 1.00000

pca_tot_df=as.data.frame(pca_tot_peak$rotation) %>% rownames_to_column(var="lib") %>% mutate(line=substr(lib,2,6))
pca_tot_df$line=as.integer(pca_tot_df$line)

map_stats=read.csv("../data/comb_map_stats_39ind.csv", header=T)

map_stat_total=map_stats %>% filter(fraction=="total")
map_stat_total$batch=as.factor(map_stat_total$batch)

pca_tot_df=pca_tot_df %>% full_join(map_stat_total, by="line")

totPCA_batch=ggplot(pca_tot_df, aes(x=PC1, y=PC2, col=batch )) + geom_point() + labs(x="PC1:0.89", y="PC2:0.043", title="Raw PAS qunatification data Total \n colored by batch ")
ggsave("../output/plots/QC_plots/TotalPCA_colBatch.png",totPCA_batch)

Saving 7 x 5 in image

totPCA_mapped=ggplot(pca_tot_df, aes(x=PC1, y=PC2, col=comb_mapped )) + geom_point() + labs(x="PC1:0.89", y="PC2:0.043", title="Raw PAS qunatification data Total \n colored by Mapped Read count")
ggsave("../output/plots/QC_plots/TotalPCA_colMapped.png",totPCA_mapped)

Saving 7 x 5 in image

pca_nuc_peak=prcomp(nuclear_Cov, center=T,scale=T)
summary(pca_nuc_peak)

Importance of components:
                          PC1     PC2     PC3     PC4     PC5     PC6
Standard deviation     5.3861 1.87775 1.62240 0.99268 0.92998 0.63513
Proportion of Variance 0.7438 0.09041 0.06749 0.02527 0.02218 0.01034
Cumulative Proportion  0.7438 0.83425 0.90174 0.92701 0.94919 0.95953
                           PC7    PC8    PC9    PC10    PC11    PC12
Standard deviation     0.53149 0.4674 0.4095 0.36160 0.32862 0.28960
Proportion of Variance 0.00724 0.0056 0.0043 0.00335 0.00277 0.00215
Cumulative Proportion  0.96677 0.9724 0.9767 0.98003 0.98280 0.98495
                          PC13    PC14   PC15    PC16    PC17    PC18
Standard deviation     0.26862 0.25414 0.2333 0.22825 0.20329 0.19277
Proportion of Variance 0.00185 0.00166 0.0014 0.00134 0.00106 0.00095
Cumulative Proportion  0.98680 0.98845 0.9899 0.99118 0.99224 0.99320
                          PC19    PC20    PC21    PC22    PC23    PC24
Standard deviation     0.18620 0.17247 0.16092 0.14244 0.13630 0.12741
Proportion of Variance 0.00089 0.00076 0.00066 0.00052 0.00048 0.00042
Cumulative Proportion  0.99409 0.99485 0.99551 0.99603 0.99651 0.99693
                          PC25    PC26    PC27    PC28    PC29    PC30
Standard deviation     0.12025 0.11377 0.11306 0.10563 0.10228 0.09219
Proportion of Variance 0.00037 0.00033 0.00033 0.00029 0.00027 0.00022
Cumulative Proportion  0.99730 0.99763 0.99796 0.99824 0.99851 0.99873
                          PC31    PC32    PC33    PC34    PC35    PC36
Standard deviation     0.08916 0.08768 0.08144 0.07916 0.07412 0.07253
Proportion of Variance 0.00020 0.00020 0.00017 0.00016 0.00014 0.00013
Cumulative Proportion  0.99893 0.99913 0.99930 0.99946 0.99960 0.99974
                          PC37    PC38    PC39
Standard deviation     0.06394 0.05721 0.05416
Proportion of Variance 0.00010 0.00008 0.00008
Cumulative Proportion  0.99984 0.99992 1.00000

pca_nuc_df=as.data.frame(pca_nuc_peak$rotation) %>% rownames_to_column(var="lib") %>% mutate(line=substr(lib,2,6))
pca_nuc_df$line=as.integer(pca_nuc_df$line)

map_stat_nuclear=map_stats %>% filter(fraction=="nuclear")
map_stat_nuclear$batch=as.factor(map_stat_nuclear$batch)

pca_nuc_df=pca_nuc_df %>% full_join(map_stat_nuclear, by="line")

nucPCA_batch=ggplot(pca_nuc_df, aes(x=PC1, y=PC2, col=batch )) + geom_point() + labs(x="PC1: 0.74", y="PC2: 0.09", title="Raw PAS qunatification data nuclear \n colored by batch ")
ggsave("../output/plots/QC_plots/NuclearPCA_colBatch.png",nucPCA_batch)

Saving 7 x 5 in image

nucPCA_mapped=ggplot(pca_nuc_df, aes(x=PC1, y=PC2, col=comb_mapped )) + geom_point() + labs(x="PC1: 0.74", y="PC2: 0.09", title="Raw PAS qunatification data nuclear \n colored by Mapped Read count")
ggsave("../output/plots/QC_plots/NuclearlPCA_colMapped.png",nucPCA_mapped)

Saving 7 x 5 in image

Q: Do the PAS read number recapitulate gene expression as it should?

Plot: scatter plot + fit (x-axis: gene TPM, y-axis: gene normalized PAS counts) total/nuclear separate

The TPM measurements come from the kalisto run I did on 18486.

tx2gene=read.table("../data/RNAkalisto/ncbiRefSeq.txn2gene.txt" ,header= F, sep="\t", stringsAsFactors = F)

txi.kallisto.tsv <- tximport("../data/RNAkalisto/abundance.tsv", type = "kallisto", tx2gene = tx2gene,countsFromAbundance="lengthScaledTPM" )

Note: importing `abundance.h5` is typically faster than `abundance.tsv`

reading in files with read_tsv

1 
removing duplicated transcript rows from tx2gene
transcripts missing from tx2gene: 99
summarizing abundance
summarizing counts
summarizing length

I need to get all of the peaks for 18486 and which gene they are in. Then I will take the gene average and divide by the number of mapped reads.

total_Cov_18486=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% filter(X18486_T>10) %>%  group_by(gene) %>% summarize(GeneSum=sum(X18486_T)) %>% mutate(GeneSumNorm=GeneSum/10.8)

#%>% mutate(NormGenePeakCov=GeneSum/10819437)

Join with the transcript TPM

TXN_abund=as.data.frame(txi.kallisto.tsv$abundance) %>% rownames_to_column(var="gene")
colnames(TXN_abund)=c("gene", "TPM")

TXN_NormGene=TXN_abund %>% inner_join(total_Cov_18486,by="gene")

Plot distribution of each variable seperatly first to understand distribution:

summary(TXN_abund$TPM)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
     0.00      0.02      1.03     36.87     14.21 101438.00 

ggplot(TXN_abund, aes(x=log10(TPM))) + geom_density(kernel="gaussian") + scale_x_log10()

Warning in self$trans$transform(x): NaNs produced

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 13505 rows containing non-finite values (stat_density).

summary(total_Cov_18486$GeneSumNorm)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
    1.019     7.963    22.963    68.164    56.551 17110.648 

ggplot(total_Cov_18486, aes(x=log10(GeneSumNorm))) + geom_density(kernel="gaussian")+ scale_x_log10()

Create a scatterplot:

TXN_NormGene=TXN_NormGene %>% filter(TPM>0) %>% filter(GeneSumNorm>0)
corr_18486Tot=ggplot(TXN_NormGene, aes(x=log10(TPM), y= log10(GeneSumNorm))) + geom_point() + labs(title="Total", x="log10 RNA seq TPM", y="log10 Peak count sum per gene")+ geom_smooth(aes(x=log10(TPM),y=log10(GeneSumNorm)),method = "lm") + annotate("text",x=5, y=5,label="R2=.42")+ geom_text(aes(label=gene),hjust=0, vjust=0)
       
corr_18486Tot       

summary(lm(log10(TPM)~log10(GeneSumNorm),TXN_NormGene)) 

Call:
lm(formula = log10(TPM) ~ log10(GeneSumNorm), data = TXN_NormGene)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.3165 -0.2546  0.0856  0.3960  2.8572 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.20648    0.01505  -13.72   <2e-16 ***
log10(GeneSumNorm)  0.92672    0.01011   91.63   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6755 on 11541 degrees of freedom
Multiple R-squared:  0.4211,    Adjusted R-squared:  0.4211 
F-statistic:  8397 on 1 and 11541 DF,  p-value: < 2.2e-16

Let me try this with the nuclear fraction:

nuclear_Cov_18486=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% filter(X18486_N>10) %>% group_by(gene) %>% summarize(GeneSum=sum(X18486_N)) %>% mutate(GeneSumNorm=GeneSum/11.4)

TXN_NormGene_Nuc=TXN_abund %>% inner_join(nuclear_Cov_18486,by="gene")

Create a scatterplot:

TXN_NormGene_Nuc=TXN_NormGene_Nuc %>% filter(TPM>0) %>% filter(GeneSumNorm>0)

corr_18486Nuc=ggplot(TXN_NormGene_Nuc, aes(x=log(TPM), y= log10(GeneSumNorm))) + geom_point() + geom_smooth(aes(x = log10(TPM +.001), y = log10(GeneSumNorm+.001)),method = "lm",se=T) + labs(title=" Nuclear", x="log10 RNA seq TPM", y="log10 Peak Sum per gene") + annotate("text",x=-3, y=5,label="R2=.32") +geom_text(aes(label=gene),hjust=0, vjust=0)
corr_18486Nuc

summary(lm(log10(TPM)~log10(GeneSumNorm),TXN_NormGene_Nuc)) 

Call:
lm(formula = log10(TPM) ~ log10(GeneSumNorm), data = TXN_NormGene_Nuc)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.7166 -0.3526  0.0734  0.4565  3.2783 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -0.13635    0.01665  -8.187 2.95e-16 ***
log10(GeneSumNorm)  0.82341    0.01102  74.727  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7629 on 11972 degrees of freedom
Multiple R-squared:  0.3181,    Adjusted R-squared:  0.318 
F-statistic:  5584 on 1 and 11972 DF,  p-value: < 2.2e-16

title <- ggdraw() + draw_label("Correlation between TPM and 3' Seq \nNA18486", fontface='bold')

plots=plot_grid(corr_18486Tot,corr_18486Nuc)

CorrelationPlot18486=plot_grid(title,plots, ncol=1 , rel_heights = c(.1,1))
ggsave(file="../output/plots/QC_plots/CorrelationWKalisto18486.png",CorrelationPlot18486)

Saving 7 x 5 in image

CorrelationPlot18486

These do not look good. We expect higher correlations. I can make variants of these plots to diagnose the problem.

topPeakCov_total_18486=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% group_by(gene) %>%  arrange(gene,desc(X18486_T)) %>% top_n(1)

Selecting by X18486_T

TXN_TopPeak=TXN_abund %>% inner_join(topPeakCov_total_18486,by="gene")


ggplot(TXN_TopPeak, aes(x=log(TPM+.001), y= log10(X18486_T+.001))) + geom_point() + geom_smooth(aes(x = log10(TPM +.001), y = log10(X18486_T+.001)),method = "lm",se=T) 

summary(lm(log10(TPM +.001)~log10(X18486_T+ .001),TXN_TopPeak)) 

Call:
lm(formula = log10(TPM + 0.001) ~ log10(X18486_T + 0.001), data = TXN_TopPeak)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5665 -0.3715  0.3014  0.6793  2.9706 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.082190   0.008386   9.801   <2e-16 ***
log10(X18486_T + 0.001) 0.360467   0.003491 103.257   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.152 on 19537 degrees of freedom
Multiple R-squared:  0.3531,    Adjusted R-squared:  0.353 
F-statistic: 1.066e+04 on 1 and 19537 DF,  p-value: < 2.2e-16

topPeakCov_nuclear_18486=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% group_by(gene) %>%  arrange(gene,desc(X18486_N)) %>% top_n(1) 

Selecting by X18486_N

TXN_TopPeak_nuc=TXN_abund %>% inner_join(topPeakCov_nuclear_18486,by="gene")


ggplot(TXN_TopPeak_nuc, aes(x=log(TPM), y= log10(X18486_N))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_N)),method = "lm",se=T) +xlim(-1,10)

Warning: Removed 3349 rows containing non-finite values (stat_smooth).

Warning: Removed 2733 rows containing missing values (geom_point).

summary(lm(log10(TPM+.001)~log10(X18486_N+.001),TXN_TopPeak_nuc)) 

Call:
lm(formula = log10(TPM + 0.001) ~ log10(X18486_N + 0.001), data = TXN_TopPeak_nuc)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4195 -0.4383  0.2779  0.6944  3.6667 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)             -0.385540   0.012683  -30.40   <2e-16 ***
log10(X18486_N + 0.001)  0.546958   0.006061   90.25   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.177 on 16302 degrees of freedom
Multiple R-squared:  0.3332,    Adjusted R-squared:  0.3331 
F-statistic:  8144 on 1 and 16302 DF,  p-value: < 2.2e-16

TXN_TopPeak_filt=TXN_TopPeak %>% filter(TPM>0) %>% filter(X18486_T>0)

ggplot(TXN_TopPeak_filt, aes(x=log(TPM), y= log10(X18486_T))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_T)),method = "lm",se=T)+ geom_text(aes(label=gene),hjust=0, vjust=0)

summary(lm(log10(TPM)~log10(X18486_T),TXN_TopPeak_filt)) 

Call:
lm(formula = log10(TPM) ~ log10(X18486_T), data = TXN_TopPeak_filt)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.3583 -0.2584  0.1128  0.4186  2.8988 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -0.906074   0.017369  -52.17   <2e-16 ***
log10(X18486_T)  0.910187   0.008167  111.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.705 on 12860 degrees of freedom
Multiple R-squared:  0.4913,    Adjusted R-squared:  0.4912 
F-statistic: 1.242e+04 on 1 and 12860 DF,  p-value: < 2.2e-16

TXN_TopPeak_filt_nuc=TXN_TopPeak_nuc %>% filter(TPM>0) %>% filter(X18486_N>0)

ggplot(TXN_TopPeak_filt_nuc, aes(x=log(TPM), y= log10(X18486_N))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_N)),method = "lm",se=T)+ geom_text(aes(label=gene),hjust=0, vjust=0)

summary(lm(log10(TPM)~log10(X18486_N),TXN_TopPeak_filt_nuc)) 

Call:
lm(formula = log10(TPM) ~ log10(X18486_N), data = TXN_TopPeak_filt_nuc)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0537 -0.3514  0.0857  0.4783  3.3509 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.021815   0.018052   -56.6   <2e-16 ***
log10(X18486_N)  0.957513   0.008878   107.8   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7692 on 13958 degrees of freedom
Multiple R-squared:  0.4545,    Adjusted R-squared:  0.4545 
F-statistic: 1.163e+04 on 1 and 13958 DF,  p-value: < 2.2e-16

outlier=c("POTEJ", "SPATA31A1", "TP53TG3B")

TXN_TopPeak_filt2_nuc= TXN_TopPeak_filt_nuc %>% filter(!(gene %in% outlier))

TXN_TopPeak_filt2= TXN_TopPeak_filt %>% filter(!(gene %in% outlier))

#nuclear 
ggplot(TXN_TopPeak_filt2_nuc, aes(x=log(TPM), y= log10(X18486_N))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_N)),method = "lm",se=T) 

summary(lm(log10(TPM)~log10(X18486_N),TXN_TopPeak_filt2_nuc)) 

Call:
lm(formula = log10(TPM) ~ log10(X18486_N), data = TXN_TopPeak_filt2_nuc)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2929 -0.3507  0.0844  0.4763  3.3486 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.01719    0.01781   -57.1   <2e-16 ***
log10(X18486_N)  0.95605    0.00876   109.1   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7588 on 13955 degrees of freedom
Multiple R-squared:  0.4605,    Adjusted R-squared:  0.4604 
F-statistic: 1.191e+04 on 1 and 13955 DF,  p-value: < 2.2e-16

#total
ggplot(TXN_TopPeak_filt2, aes(x=log(TPM), y= log10(X18486_T))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_T)),method = "lm",se=T)

summary(lm(log10(TPM)~log10(X18486_T),TXN_TopPeak_filt)) 

Call:
lm(formula = log10(TPM) ~ log10(X18486_T), data = TXN_TopPeak_filt)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.3583 -0.2584  0.1128  0.4186  2.8988 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -0.906074   0.017369  -52.17   <2e-16 ***
log10(X18486_T)  0.910187   0.008167  111.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.705 on 12860 degrees of freedom
Multiple R-squared:  0.4913,    Adjusted R-squared:  0.4912 
F-statistic: 1.242e+04 on 1 and 12860 DF,  p-value: < 2.2e-16

OnePeak_nuc=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% group_by(gene) %>% tally() %>% filter(n==1) %>% select(gene)

OnePeak_nuc_cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% inner_join(OnePeak_nuc, by="gene") %>% inner_join(TXN_abund, by="gene")

OnePeak_nuc_cov=OnePeak_nuc_cov %>% filter(TPM>0) %>% filter(X18486_N>0)

ggplot(OnePeak_nuc_cov, aes(x=log10(TPM), y= log10(X18486_N))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_N)),method = "lm",se=T) + geom_text(aes(label=gene),hjust=0, vjust=0)

summary(lm(log10(TPM)~log10(X18486_N),OnePeak_nuc_cov)) 

Call:
lm(formula = log10(TPM) ~ log10(X18486_N), data = OnePeak_nuc_cov)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.8689 -0.5502  0.0730  0.5916  3.2499 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.17296    0.05663  -20.71   <2e-16 ***
log10(X18486_N)  0.84559    0.04475   18.90   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9359 on 801 degrees of freedom
Multiple R-squared:  0.3084,    Adjusted R-squared:  0.3075 
F-statistic: 357.1 on 1 and 801 DF,  p-value: < 2.2e-16

OnePeakN_outlier=c("POTEJ")


OnePeak_nuc_cov_filt= OnePeak_nuc_cov %>% filter(!(gene %in% OnePeakN_outlier))


ggplot(OnePeak_nuc_cov_filt, aes(x=log10(TPM), y= log10(X18486_N))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_N)),method = "lm",se=T) + geom_text(aes(label=gene),hjust=0, vjust=0)

summary(lm(log10(TPM)~log10(X18486_N),OnePeak_nuc_cov_filt)) 

Call:
lm(formula = log10(TPM) ~ log10(X18486_N), data = OnePeak_nuc_cov_filt)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2074 -0.5615  0.0608  0.5899  3.2384 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.14670    0.05345  -21.45   <2e-16 ***
log10(X18486_N)  0.83081    0.04221   19.68   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8823 on 800 degrees of freedom
Multiple R-squared:  0.3263,    Adjusted R-squared:  0.3254 
F-statistic: 387.4 on 1 and 800 DF,  p-value: < 2.2e-16

OnePeak_tot=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% group_by(gene) %>% tally() %>% filter(n==1) %>% select(gene)

OnePeak_tot_cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% inner_join(OnePeak_tot, by="gene") %>% inner_join(TXN_abund, by="gene")

OnePeak_tot_cov=OnePeak_tot_cov %>% filter(TPM>0) %>% filter(X18486_T>0)

ggplot(OnePeak_tot_cov, aes(x=log10(TPM), y= log10(X18486_T))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_T)),method = "lm",se=T) + geom_text(aes(label=gene),hjust=0, vjust=0)

summary(lm(log10(TPM)~log10(X18486_T),OnePeak_tot_cov)) 

Call:
lm(formula = log10(TPM) ~ log10(X18486_T), data = OnePeak_tot_cov)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.5586 -0.5457  0.1636  0.6173  2.5539 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.22867    0.07551  -16.27   <2e-16 ***
log10(X18486_T)  0.86571    0.04861   17.81   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9762 on 536 degrees of freedom
Multiple R-squared:  0.3718,    Adjusted R-squared:  0.3706 
F-statistic: 317.2 on 1 and 536 DF,  p-value: < 2.2e-16

OnePeak_tot_cov_filt= OnePeak_tot_cov %>% filter(!(gene %in% OnePeakN_outlier))


ggplot(OnePeak_tot_cov_filt, aes(x=log10(TPM), y= log10(X18486_T))) + geom_point() + geom_smooth(aes(x = log10(TPM), y = log10(X18486_T)),method = "lm",se=T) + geom_text(aes(label=gene),hjust=0, vjust=0)

summary(lm(log10(TPM)~log10(X18486_T),OnePeak_tot_cov_filt)) 

Call:
lm(formula = log10(TPM) ~ log10(X18486_T), data = OnePeak_tot_cov_filt)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4792 -0.5680  0.1435  0.6233  2.5532 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.17715    0.07013  -16.79   <2e-16 ***
log10(X18486_T)  0.83818    0.04510   18.59   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9039 on 535 degrees of freedom
Multiple R-squared:  0.3923,    Adjusted R-squared:  0.3912 
F-statistic: 345.4 on 1 and 535 DF,  p-value: < 2.2e-16

Q: For each gene, what percentage of reads assigned fall within 1, 2, 3, etc… peaks, we would expect that for many genes >90% of the reads fall within 1 peak, for a few 2 peaks, etc…?

Plot: Y-axis: Number of genes, X-axis: how many peaks is needed to “capture” 90%, 80%, … 50% of the reads assigned to that gene (using different colors).

Start with analysis to see how many peaks are needed to capture 90% of the reads assigned to the gene. I will start by looking at the number of reads that map to peaks in genes. To do this I can group on genes in the peak coverage and get the sum.

nuclear_covBygene=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% group_by(gene) %>% summarize(GeneSum=sum(X18486_N)) %>% mutate(per90=GeneSum*.9)%>% mutate(per80=GeneSum*.8)%>% mutate(per70=GeneSum*.7)%>% mutate(per60=GeneSum*.6)%>% mutate(per50=GeneSum*.5) 
total_covBygene=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% group_by(gene) %>% summarize(GeneSum=sum(X18486_T))%>% mutate(per90=GeneSum*.9)%>% mutate(per80=GeneSum*.8)%>% mutate(per70=GeneSum*.7)%>% mutate(per60=GeneSum*.6)%>% mutate(per50=GeneSum*.5)

Write these out to use them in the script:

write.table(file="../data/UnderstandPeaksQC/Nuclear_PerCovbyGene.txt", nuclear_covBygene, quote=F, col.names = T, row.names = F)
write.table(file="../data/UnderstandPeaksQC/Total_PerCovbyGene.txt", total_covBygene, quote=F, col.names = T,row.names = F)

Try in R

#groupedNuclear=sumforeachgene %>% sort(ind) %>% cumulativesum %>% dividebygenesum %>% filter(only90) %>% count()
#remove genes with 0 count sum 

nuclear_90Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% group_by(gene) %>%  arrange(gene,desc(X18486_N)) %>%  mutate(SUM = cumsum(X18486_N)) %>% full_join(nuclear_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.9) %>% tally() %>% rename("90"=n)

nuclear_80Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% group_by(gene) %>%  arrange(gene,desc(X18486_N)) %>%  mutate(SUM = cumsum(X18486_N)) %>% full_join(nuclear_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.8) %>% tally() %>% rename("80"=n)

nuclear_70Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% group_by(gene) %>%  arrange(gene,desc(X18486_N)) %>%  mutate(SUM = cumsum(X18486_N)) %>% full_join(nuclear_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.7) %>% tally() %>% rename("70"=n)

nuclear_60Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% group_by(gene) %>%  arrange(gene,desc(X18486_N)) %>%  mutate(SUM = cumsum(X18486_N)) %>% full_join(nuclear_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.6) %>% tally() %>% rename("60"=n)

nuclear_50Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Nuclear_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_N) %>% group_by(gene) %>%  arrange(gene,desc(X18486_N)) %>%  mutate(SUM = cumsum(X18486_N)) %>% full_join(nuclear_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.5) %>% tally() %>% rename("50"=n)

Join these to plot them:

nuclear_PercentPeakCov= nuclear_90Cov %>% left_join(nuclear_80Cov, by="gene") %>% left_join(nuclear_70Cov, by="gene") %>% left_join(nuclear_60Cov, by="gene") %>% left_join(nuclear_50Cov, by="gene")

nuclear_PercentPeakCov_melt=melt(nuclear_PercentPeakCov,id.vars = "gene")

nucPeakCov=ggplot(nuclear_PercentPeakCov_melt, aes(x=value,fill=variable))+ geom_histogram(position="dodge", bins=30) + labs(y="Number of Genes", x="Number of Peaks", title="Nuclear: Number of Peaks to capture % of Gene count") + facet_grid(~variable) + xlim(0,30)


nucPeakCov_cdf=ggplot(nuclear_PercentPeakCov_melt, aes(x=value,col=variable))+ stat_ecdf(geom="step")+ labs(y="Percent of Genes", x="Number of Peaks", title="Nuclear: Number of Peaks to capture % of Gene count") + scale_x_continuous(breaks=seq(1,30,2),limits=c(0,30))

ggplot(nuclear_PercentPeakCov_melt, aes(x=value,fill=variable, by=variable))+ geom_density(alpha=.4) + labs(y="Number of Genes", x="Number of Peaks", title="Nuclear: Number of Peaks to capture % of Gene count")

Try this with total:

total_90Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% group_by(gene) %>%  arrange(gene,desc(X18486_T)) %>%  mutate(SUM = cumsum(X18486_T)) %>% full_join(total_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.9) %>% tally() %>% rename("90"=n)

total_80Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% group_by(gene) %>%  arrange(gene,desc(X18486_T)) %>%  mutate(SUM = cumsum(X18486_T)) %>% full_join(total_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.8) %>% tally() %>% rename("80"=n)

total_70Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% group_by(gene) %>%  arrange(gene,desc(X18486_T)) %>%  mutate(SUM = cumsum(X18486_T)) %>% full_join(total_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.7) %>% tally() %>% rename("70"=n)

total_60Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% group_by(gene) %>%  arrange(gene,desc(X18486_T)) %>%  mutate(SUM = cumsum(X18486_T)) %>% full_join(total_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.6) %>% tally() %>% rename("60"=n)

total_50Cov=read.table("../data/PeakCounts/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.Total_fixed.fc", header=T, stringsAsFactors = F)[,1:7] %>% separate(Geneid, into=c("peak", "chr", "start", "end", "strand", "gene"), sep=":") %>% select(gene, X18486_T) %>% group_by(gene) %>%  arrange(gene,desc(X18486_T)) %>%  mutate(SUM = cumsum(X18486_T)) %>% full_join(total_covBygene,by="gene") %>% filter(GeneSum >0) %>% mutate(perSum=SUM/GeneSum) %>% mutate(perSum_lag=lag(perSum,1)) %>%  replace_na(list(perSum_lag =0)) %>% filter(perSum_lag<.5) %>% tally() %>% rename("50"=n)

Put together:

total_PercentPeakCov= total_90Cov %>% left_join(total_80Cov, by="gene") %>% left_join(total_70Cov, by="gene") %>% left_join(total_60Cov, by="gene") %>% left_join(total_50Cov, by="gene")

total_PercentPeakCov_melt=melt(total_PercentPeakCov,id.vars = "gene")

totPeakCov=ggplot(total_PercentPeakCov_melt, aes(x=value,fill=variable))+ geom_histogram(position="dodge", bins=30) + labs(y="Number of Genes", x="Number of Peaks", title="Total: Number of Peaks to capture % of Gene count") + facet_grid(~variable) + xlim(0,30)


totPeakCov_cdf=ggplot(total_PercentPeakCov_melt, aes(x=value,col=variable))+ stat_ecdf(geom="step")+ labs(y="Percent of Genes", x="Number of Peaks", title="Total: Number of Peaks to capture % of Gene count") + scale_x_continuous(breaks=seq(1,30,2),limits=c(0,30))




ggplot(total_PercentPeakCov_melt, aes(x=value,fill=variable, by=variable))+ geom_density(alpha=.4) + labs(y="Number of Genes", x="Number of Peaks", title="Nuclear: Number of Peaks to capture % of Gene count")

PeakCovPerGeneCount=plot_grid(totPeakCov,nucPeakCov, ncol = 1)

Warning: Removed 85 rows containing non-finite values (stat_bin).

Warning: Removed 810 rows containing non-finite values (stat_bin).

ggsave(file="../output/plots/QC_plots/PeakCovPerGeneCount.png",PeakCovPerGeneCount)

Saving 7 x 5 in image

PeakCovPerGeneCountCDF=plot_grid(totPeakCov_cdf,nucPeakCov_cdf,ncol=1)

Warning: Removed 85 rows containing non-finite values (stat_ecdf).

Warning: Removed 810 rows containing non-finite values (stat_ecdf).

ggsave(file="../output/plots/QC_plots/PeakCovPerGeneCountCDF.png",PeakCovPerGeneCountCDF)

Saving 7 x 5 in image

These look good. We could use this to filter. I would compute the percent of the gene each peak covers in each individual then I could filter peaks that cover % t() fc_peaks=as.data.frame(fc_peaks) colnames(fc_peaks)=as.character(unlist(fc_peaks[1,])) fc_peaks=fc_peaks[-1,] fc_peaks\(Assigned=as.numeric(as.character(fc_peaks\)Assigned)) fc_peaks\(Unassigned_NoFeatures=as.numeric(as.character(fc_peaks\)Unassigned_NoFeatures)) ```

I need to separate the libraries by line and fraction.

fc_peaks=fc_peaks %>% separate(Status, into=c("line", "fraction"), sep="_") %>% mutate(PerReadPeak=Assigned/(Assigned+Unassigned_NoFeatures))

This number is the reads assigned to peaks out of all reads mapping to genome.

I can now melt these data by line and fraction

fc_peaks_melt=melt(fc_peaks, id.vars = c("line", "fraction"))

Warning: attributes are not identical across measure variables; they will
be dropped

fc_peaks_melt_PerRead=fc_peaks_melt %>% filter(variable=="PerReadPeak")
fc_peaks_melt_PerRead$value=as.numeric(fc_peaks_melt_PerRead$value)

ggplot(fc_peaks_melt_PerRead,aes( x=line, y=value, by=fraction, fill=fraction))+ geom_col(pos="dodge") +theme(axis.text.x = element_text(angle = 90, hjust = 1),axis.text.y = element_text(size=12),axis.title.y=element_text(size=10,face="bold"), axis.title.x=element_text(size=12,face="bold"))+ scale_fill_manual(values=c("deepskyblue3","darkviolet")) + labs(title="Percent of reads mapping to peaks by line and fraction", y="Reads mapping to peaks/all mapping reads")

It may be more interesting to look at this by fraction, with error bars.

fc_peaks_melt_PerRead_byfrac= fc_peaks_melt_PerRead %>% group_by(fraction) %>% summarise(mean=mean(value), sd=sd(value))

Plot this:

ggplot(fc_peaks_melt_PerRead_byfrac,aes(x=fraction, y=mean, fill=fraction)) + geom_col()+ geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)+ theme(axis.text.y = element_text(size=12),axis.title.y=element_text(size=10,face="bold"), axis.title.x=element_text(size=12,face="bold"))+ scale_fill_manual(values=c("deepskyblue3","darkviolet"))+ labs(title="Percent of reads mapping to peaks by fraction", y="Reads mapping to peaks/all mapping reads")

Now I want to look at how many reads map to gene. I will use the transcript annotations that I used for the peaks.

• /project2/gilad/briana/genome_anotation_data/ncbiRefSeq_sm_noChr.sort.mRNA.bed

I need to make this an SAF file.
* GeneID * Chr * Start * End * Strand

RefSeqmRNA2SAF.py

#python
from misc_helper import *
fout = file("/project2/gilad/briana/genome_anotation_data/ncbiRefSeq_sm_noChr.sort.mRNA.SAF","w")
fout.write("GeneID\tChr\tStart\tEnd\tStrand\n")
for ln in open("/project2/gilad/briana/genome_anotation_data/ncbiRefSeq_sm_noChr.sort.mRNA.bed"):
    chrom, start, end, gene, score, strand = ln.split()
    start_i=int(start)
    end_i=int(end)
    fout.write("%s\t%s\t%d\t%d\t%s\n"%(gene, chrom, start_i, end_i, strand))
fout.close()

ref_geneTranscript_fc.sh

#!/bin/bash

#SBATCH --job-name=ref_geneTranscript_fc
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=ref_geneTranscript_fc.out
#SBATCH --error=ref_geneTranscript_fc.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END


module load Anaconda3
source activate three-prime-env

featureCounts -O -a /project2/gilad/briana/genome_anotation_data/ncbiRefSeq_sm_noChr.sort.mRNA.SAF -F SAF -o /project2/gilad/briana/threeprimeseq/data/UnderstandPeaksQC/RefSeqTranscript_AllLibraries.fc /project2/gilad/briana/threeprimeseq/data/sort/*sort.bam -s 2

fix_Genefc_summary.py

infile= open("/Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/UnderstandPeaksQC/RefSeqTranscript_AllLibraries.fc.summary", "r")
fout = open("/Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/UnderstandPeaksQC/RefSeqTranscript_AllLibraries.fc.summary_fixed",'w')
for line, i in enumerate(infile):
    if line == 0:
        i_list=i.split()
        libraries=[i_list[0]]
        for sample in i_list[1:]:
            full = sample.split("/")[7]
            samp= full.split("-")[2:4]
            lim="_"
            samp_st=lim.join(samp)
            libraries.append(samp_st)
        print(libraries)
        first_line= "\t".join(libraries)
        fout.write(first_line + '\n' )
    else:
        fout.write(i)
fout.close()

fc_gene_peaks=read.table("../data/UnderstandPeaksQC/RefSeqTranscript_AllLibraries.fc.summary_fixed", stringsAsFactors = F) %>% t()
fc_gene_peaks=as.data.frame(fc_gene_peaks)
colnames(fc_gene_peaks)=as.character(unlist(fc_gene_peaks[1,]))
fc_gene_peaks=fc_gene_peaks[-1,]
fc_gene_peaks$Assigned=as.numeric(as.character(fc_gene_peaks$Assigned))
fc_gene_peaks$Unassigned_NoFeatures=as.numeric(as.character(fc_gene_peaks$Unassigned_NoFeatures))

I need to separate the libraries by line and fraction.

fc_gene_peaks=fc_gene_peaks %>% separate(Status, into=c("line", "fraction"), sep="_") %>% mutate(PerReadPeak=Assigned/(Assigned+Unassigned_NoFeatures))

Melt this:

fc_gene_peaks_melt=melt(fc_gene_peaks, id.vars = c("line", "fraction"))

Warning: attributes are not identical across measure variables; they will
be dropped

fc_gene_peaks_PerRead=fc_gene_peaks_melt %>% filter(variable=="PerReadPeak")
fc_gene_peaks_PerRead$value=as.numeric(fc_gene_peaks_PerRead$value)

GGplot:

ggplot(fc_gene_peaks_PerRead,aes( x=line, y=value, by=fraction, fill=fraction))+ geom_col(pos="dodge") +theme(axis.text.x = element_text(angle = 90, hjust = 1),axis.text.y = element_text(size=12),axis.title.y=element_text(size=10,face="bold"), axis.title.x=element_text(size=12,face="bold"))+ scale_fill_manual(values=c("deepskyblue3","darkviolet")) + labs(title="Percent of reads mapping to Transcripts by line and fraction", y="Reads mapping to transcripts/all mapping reads")

Do this by fraction.

fc_gene_peaks_PerRead_byfrac= fc_gene_peaks_PerRead %>% group_by(fraction) %>% summarise(mean=mean(value), sd=sd(value))

Plot this:

ggplot(fc_gene_peaks_PerRead_byfrac,aes(x=fraction, y=mean, fill=fraction)) + geom_col()+ geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)+ theme(axis.text.y = element_text(size=12),axis.title.y=element_text(size=10,face="bold"), axis.title.x=element_text(size=12,face="bold"))+ scale_fill_manual(values=c("deepskyblue3","darkviolet"))+ labs(title="Percent of reads mapping to Transcripts by fraction", y="Reads mapping to Transcripts/all mapping reads")

It would be nice to have this in one plot. In order to do this I want to join the PerReadPeak from both and melt. this way the variable can be peak or transcript.

fc_peaks_sel=fc_peaks %>% select(c("line", "fraction", "PerReadPeak"))

fc_gene_peaks_sel=fc_gene_peaks %>% select(c("line", "fraction", "PerReadPeak"))

fcGene_and_Transcript=fc_peaks_sel %>% left_join(fc_gene_peaks_sel, by=c("line","fraction"))

colnames(fcGene_and_Transcript)=c("Line", "Fraction", "Peaks", "Genes")


fcGene_and_Transcript_melt=melt(fcGene_and_Transcript, id.vars=c("Line","Fraction"))


fcGene_and_Transcript_melt_sum=fcGene_and_Transcript_melt %>% group_by(Fraction,variable) %>% summarise(mean=mean(value), sd=sd(value))

reads2featuresPlot=ggplot(fcGene_and_Transcript_melt_sum,aes(x=Fraction, y=mean, fill=Fraction)) + geom_col()+ geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)+ theme(axis.text.y = element_text(size=12),axis.title.y=element_text(size=10,face="bold"), axis.title.x=element_text(size=12,face="bold"))+ scale_fill_manual(values=c("deepskyblue3","darkviolet"))+ labs(title="Percent of reads mapping to feature by fraction", y="Reads mapping to Feature/all mapping reads") + facet_grid(~variable)

reads2featuresPlot

ggsave(file="../output/plots/QC_plots/reads2featuresPlot.png", reads2featuresPlot)

Saving 7 x 5 in image

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  tximport_1.8.0  devtools_1.13.6 reshape2_1.4.3 
 [5] cowplot_0.9.3   workflowr_1.1.1 forcats_0.3.0   stringr_1.3.1  
 [9] dplyr_0.7.6     purrr_0.2.5     readr_1.1.1     tidyr_0.8.1    
[13] tibble_1.4.2    ggplot2_3.0.0   tidyverse_1.2.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] memoise_1.1.0     evaluate_0.11     labeling_0.3     
[25] knitr_1.20        broom_0.5.0       Rcpp_0.12.19     
[28] scales_1.0.0      backports_1.1.2   jsonlite_1.5     
[31] hms_0.4.2         digest_0.6.17     stringi_1.2.4    
[34] grid_3.5.1        rprojroot_1.3-2   cli_1.0.1        
[37] tools_3.5.1       magrittr_1.5      lazyeval_0.2.1   
[40] crayon_1.3.4      whisker_0.3-2     pkgconfig_2.0.2  
[43] xml2_1.2.0        lubridate_1.7.4   assertthat_0.2.0 
[46] rmarkdown_1.10    httr_1.3.1        rstudioapi_0.8   
[49] R6_2.3.0          nlme_3.1-137      git2r_0.23.0     
[52] compiler_3.5.1