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(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=TPM)) + geom_density(kernel="gaussian") + scale_x_log10()

Warning: Transformation introduced infinite values in continuous x-axis

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

summary(total_Cov_18486$GeneSum)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
    11.0     86.0    248.0    736.2    610.8 184795.0 

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

Create a scatterplot:

corr_18486Tot=ggplot(TXN_NormGene, aes(x=log10(TPM+ .001), y= log10(GeneSum+.001))) + geom_point() + labs(title="Total", x="log10 RNA seq TPM", y="log10 Peak count sum per gene")+ geom_smooth(aes(x=log10(TPM +.001),y=log10(GeneSum+ .001)),method = "lm") + xlim(-2,6) + annotate("text",x=5, y=5,label="R2=.36")
       
corr_18486Tot       

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

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

summary(lm(log10(TPM +.001)~log10(GeneSum+ .001),TXN_NormGene)) 

Call:
lm(formula = log10(TPM + 0.001) ~ log10(GeneSum + 0.001), data = TXN_NormGene)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.8203 -0.2252  0.1609  0.5026  3.1060 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -1.66435    0.03240  -51.37   <2e-16 ***
log10(GeneSum + 0.001)  1.08820    0.01324   82.17   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9111 on 11954 degrees of freedom
Multiple R-squared:  0.3609,    Adjusted R-squared:  0.3609 
F-statistic:  6752 on 1 and 11954 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(NormGenePeakCov=GeneSum/11405271)

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

Create a scatterplot:

corr_18486Nuc=ggplot(TXN_NormGene_Nuc, aes(x=log(TPM+.001), y= log10(GeneSum+.001))) + geom_point() + geom_smooth(aes(x = log10(TPM +.001), y = log10(GeneSum+.001)),method = "lm",se=T) + labs(title=" Nuclear", x="log10 RNA seq TPM", y="log10 Peak Sum per gene") +xlim(-5,10) + annotate("text",x=-3, y=5,label="R2=.30")
corr_18486Nuc

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

summary(lm(log10(TPM +.001)~log10(GeneSum+ .001),TXN_NormGene_Nuc)) 

Call:
lm(formula = log10(TPM + 0.001) ~ log10(GeneSum + 0.001), data = TXN_NormGene_Nuc)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.3659 -0.3288  0.1401  0.5747  3.5594 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -1.58670    0.03432  -46.23   <2e-16 ***
log10(GeneSum + 0.001)  1.01005    0.01380   73.21   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9916 on 12455 degrees of freedom
Multiple R-squared:  0.3009,    Adjusted R-squared:  0.3008 
F-statistic:  5360 on 1 and 12455 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)

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

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

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

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

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)

Python method
I think the next step will be best accomplished with a python script where I can go line by line. I want to order the peak counts in a dictionary. This will be double dictionary because I need gene and individual. I want the key to be the gene and the value to be an ordered list of the peak counts. I can ask how many elements it took to acieve the % read coverage.

1. Build dictionary

2. Refer to the *_covBygene gene number pairs to ask how many peaks get you to that number.

3. Write out the gene and the number to capture the % of interest

4. The file should take a percentage between 50 and 90 and total/nuclear. I can then run the script on each value and join the dataframes in R.

5. filter genes with 0 reads

File is code/PeaksToCoverPerReads.py. I am stuck with the dictionary creation because the same values are there for each gene. All values are being added. The gene key is not being recognized.

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)

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)


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

Q: What % of reads are assigned to a peak? Of these, what % of reads are assigned to a gene?

Within 50bp of an exon (more relevant for total)?

I want to know the percent of reads that are assigned to our peaks. I can get this information from the peak feature counts summaries. In order to look at the reads assigned to genes I will need to use feature counts with the gene annotation file.

• File: /project2/gilad/briana/threeprimeseq/data/filtPeakOppstrand_cov/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.fc.summary

Feature count takes in the bam files and an SAF annotation. For this one I used the peaks woth the transcript level annotation. I will fix the column names with python.

fix_fc_summary.py

infile= open("/Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/UnderstandPeaksQC/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.fc.summary", "r")
fout = open("/Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/UnderstandPeaksQC/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.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()

I care about Unassigned_NoFeatures and Assigned. These numbers add to the number of reads that map to the genome.

fc_peaks=read.table("../data/UnderstandPeaksQC/filtered_APApeaks_merged_allchrom_refseqGenes.Transcript_sm_quant.fc.summary_fixed", stringsAsFactors = F) %>% 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