PacPio Reads in my peaks

Look at published pacbio

library(workflowr)

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

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.4.0
✔ readr   1.1.1     ✔ forcats 0.3.0

Warning: package 'stringr' was built under R version 3.5.2

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

library(reshape2)

Attaching package: 'reshape2'

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

    smiths

library(cowplot)

Attaching package: 'cowplot'

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

    ggsave

Ankeeta has been working with 3 pac bio libraries for whole LCLs. The meged bam file has 4,164,259 reads. I want to look at how many of these reads cover my peaks. It would be best to know how many reads ends

I need to fix the strand for my peaks and give them to her.

fixPeaks4Ankeeta.py

In=open("/project2/gilad/briana/threeprimeseq/data/mergedPeaks_noMP_GeneLoc/Filtered_APApeaks_merged_allchrom_noMP.sort.named.noCHR_geneLocParsed_withAnno.SAF","r")
Out="/project2/yangili1/PAPeaks_STARMap_GeneLocAnno.bed"

def fix_strand(Fin,Fout):
    fout=open(Fout,"w")
    for n, ln in enumerate(Fin):
        if n == 0: 
            continue
        else: 
            id, chrom, start, end, strand = ln.split()
            if strand=="+":
                chromF="chr" + chrom
                peak=id.split(":")[0]
                geneLoc=id.split(":")[5:]
                geneLocF=":".join(geneLoc)
                newID=peak + ":" + geneLocF
                score="."
                fout.write("%s\t%s\t%s\t%s\t%s\t-\n"%(chromF,start,end,newID,score))
            else:
                chromF="chr" + chrom
                peak=id.split(":")[0]
                geneLoc=id.split(":")[5:]
                geneLocF=":".join(geneLoc)
                newID=peak + ":" + geneLocF
                score="."
                fout.write("%s\t%s\t%s\t%s\t%s\t+\n"%(chromF,start,end,newID,score))
    fout.close()
    
    
fix_strand(In, Out)

Add average usage to this:

use similar code to filter_5percUsagePeaks.R

counts only numeric are in /project2/gilad/briana/threeprimeseq/data/phenotypes_filtPeakTranscript_noMP_GeneLocAnno/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Total.fixed.pheno.CountsOnlyNumeric.txt I will take the mean for each row of this and use it as the score in the bed file.

Run this interactively

library(dplyr)
totUsage=read.table("/project2/gilad/briana/threeprimeseq/data/phenotypes_filtPeakTranscript_noMP_GeneLocAnno/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Total.fixed.pheno.CountsOnlyNumeric.txt", header=F)
peakBed=read.table("/project2/yangili1/PAPeaks_STARMap_GeneLocAnno.bed", header=F, col.names = c("chr", "start", "end", "ID", "score", "strand"), stringsAsFactors = F)


MeanUsage=rowMeans(totUsage)

outBed=as.data.frame(cbind(peakBed, MeanUsage)) %>% select(chr, start, end, ID, MeanUsage, strand)

write.table(outBed,file="/project2/yangili1/PAPeaks_STARMap_GeneLocAnno_withMeanUsage.bed", row.names=F, col.names=F, quote = F, sep="\t")  

Result from pac bio overlap:

• /project2/yangili1/ankeetashah/APA_tools/peaks/threeprime_noMP_pacbio_UTR.coverage
• project2/yangili1/ankeetashah/APA_tools/peaks/threeprime_noMP_pacbio_intron.coverage

Make some plots for this: Distribution of reads ending at each peak

covNames=c("chr", "start", "end", "ID", "score", "strand", "cov")
utrCov=read.table("../data/pacbio_cov/threeprime_noMP_pacbio_UTR.coverage", stringsAsFactors = F, col.names = covNames)
intronCov=read.table("../data/pacbio_cov/threeprime_noMP_pacbio_intron.coverage", stringsAsFactors = F,col.names = covNames) %>% separate(ID, into=c("peak", "gene","loc"), sep=":") %>% filter(loc=="intron")

Plot the distributions:

ggplot(utrCov, aes(x=log10(cov + 1))) + geom_density() + labs(title="PacBio reads ending at each UTR peak", x="log10(nReads+1)")

Past versions of unnamed-chunk-5-1.png

ggplot(intronCov, aes(x=log10(cov + 1))) + geom_density() + labs(title="PacBio reads ending at each intronic peak", x="log10(nReads+1)")

Past versions of unnamed-chunk-6-1.png

summary(intronCov$cov)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
   0.000    0.000    0.000    0.611    0.000 4145.000 

summary(utrCov$cov)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
    0.00     1.00     5.00    24.93    18.00 15301.00 

Proportion of peaks with coverage

intronCov_0=intronCov %>% filter(cov > 0) %>% nrow()/ nrow(intronCov)
intronCov_10=intronCov %>% filter(cov >= 10) %>% nrow()/nrow(intronCov)
intronCov_100=intronCov %>% filter(cov >= 100) %>% nrow()/nrow(intronCov)
intronCov_1000=intronCov %>% filter(cov >= 1000) %>% nrow()/nrow(intronCov)



utrCov_0=utrCov %>% filter(cov > 0) %>% nrow() / nrow(utrCov)
utrCov_10=utrCov %>% filter(cov >= 10) %>% nrow()/ nrow(utrCov)
utrCov_100=utrCov %>% filter(cov >= 100) %>% nrow()/ nrow(utrCov)
utrCov_1000=utrCov %>% filter(cov >= 1000) %>% nrow()/ nrow(utrCov)


Reads=c("1 read", "10 reads", "100 reads", "1000 reads")
UTR=c(utrCov_0, utrCov_10, utrCov_100, utrCov_1000)
Intron=c(intronCov_0,intronCov_10,intronCov_100,intronCov_1000)

covDF=as.data.frame(cbind(Reads,UTR,Intron))
covDF$UTR=as.numeric(as.character(covDF$UTR))
covDF$Intron=as.numeric(as.character(covDF$Intron))


covDF_melt=melt(covDF, id.vars = "Reads")
colnames(covDF_melt)=c("ReadCutoff", "Location", "Proportion" )

propcovbycutff=ggplot(covDF_melt, aes(x=ReadCutoff, fill=Location, y=Proportion, by=Location)) + geom_bar(position="dodge",stat="identity" ) + labs(y="Proportion of PAS covered", x="At least X reads ending at PAS", title="Proportion of PAS with PacBio read ending in it") + scale_fill_manual(values=c("blue","red")) + annotate("text", label="UTR PAS = 29,687", x="1000 reads", y=.7) + annotate("text", label="Intronic PAS = 87,733", x="1000 reads", y=.6)

propcovbycutff

Past versions of unnamed-chunk-9-1.png

ggsave(propcovbycutff, filename = "../output/plots/PacBioReadsEndingAtPAS.png")

Saving 7 x 5 in image

Subset for peaks that passed the 5% coverage cutoff.

tot5Perc=read.table("../data/PeakUsage_noMP_GeneLocAnno/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno.NoMP_sm_quant.Total_fixed.pheno.5percPeaks.txt", stringsAsFactors = F, col.names = c("chr", "start", "end", "gene", "strand", "peak", "avgUsage"))
nuc5Perc=read.table("../data/PeakUsage_noMP_GeneLocAnno/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno.NoMP_sm_quant.Nuclear_fixed.pheno.5percPeaks.txt",col.names = c("chr", "start", "end", "gene", "strand", "peak", "avgUsage"),stringsAsFactors = F)

intronCov_tot5perc= intronCov %>% semi_join(tot5Perc, by="peak")

intronCov_nuc5perc= intronCov %>% semi_join(nuc5Perc, by="peak")

intronCov_tot5perc_0= intronCov_tot5perc %>% filter(cov> 0) %>% nrow()/ nrow(intronCov_tot5perc)
intronCov_tot5perc_10= intronCov_tot5perc %>% filter(cov>= 10) %>% nrow()/ nrow(intronCov_tot5perc)
intronCov_tot5perc_100= intronCov_tot5perc %>% filter(cov>= 100) %>% nrow()/ nrow(intronCov_tot5perc)
intronCov_tot5perc_1000= intronCov_tot5perc %>% filter(cov>= 1000) %>% nrow()/ nrow(intronCov_tot5perc)


intronCov_nuc5perc_0= intronCov_nuc5perc %>% filter(cov> 0) %>% nrow()/ nrow(intronCov_nuc5perc)
intronCov_nuc5perc_10= intronCov_nuc5perc %>% filter(cov>= 10) %>% nrow()/ nrow(intronCov_nuc5perc)
intronCov_nuc5perc_100= intronCov_nuc5perc %>% filter(cov>= 100) %>% nrow()/ nrow(intronCov_nuc5perc)
intronCov_nuc5perc_1000= intronCov_nuc5perc %>% filter(cov>= 1000) %>% nrow()/ nrow(intronCov_nuc5perc)

Reads=c("1 read", "10 reads", "100 reads", "1000 reads")
Total=c(intronCov_tot5perc_0, intronCov_tot5perc_10, intronCov_tot5perc_100, intronCov_tot5perc_1000)
Nuclear=c(intronCov_nuc5perc_0,intronCov_nuc5perc_10,intronCov_nuc5perc_100,intronCov_nuc5perc_1000)

cov_5perDF=as.data.frame(cbind(Reads,Total,Nuclear))
cov_5perDF$Total=as.numeric(as.character(cov_5perDF$Total))
cov_5perDF$Nuclear=as.numeric(as.character(cov_5perDF$Nuclear))

cov_5perDF_melt=melt(cov_5perDF, id.vars = "Reads")
colnames(cov_5perDF_melt)=c("ReadCutoff", "Fraction", "Proportion" )

ggplot(cov_5perDF_melt, aes(x=ReadCutoff, fill=Fraction, y=Proportion, by=Fraction)) + geom_bar(position="dodge",stat="identity" ) + labs(y="Proportion of PAS covered", x="At least X reads ending at PAS", title="Proportion of PAS with PacBio read ending in it \n Intron PAS with 5% mean Usage") + scale_fill_brewer(palette = "Set1") + annotate("text", label="Total PAS = 16,662", x="1000 reads", y=.28) + annotate("text", label="Nuclear PAS = 18,829", x="1000 reads", y=.25)

Past versions of unnamed-chunk-13-1.png

Do this for the UTR

UTRCov_tot5perc= utrCov %>% separate(ID, into=c("peak", "gene","loc"), sep=":") %>% semi_join(tot5Perc, by="peak")

UTRCov_nuc5perc= utrCov %>% separate(ID, into=c("peak", "gene","loc"), sep=":") %>% semi_join(nuc5Perc, by="peak")

utrCov_tot5perc_0= UTRCov_tot5perc %>% filter(cov> 0) %>% nrow()/ nrow(UTRCov_tot5perc)
utrCov_tot5perc_10= UTRCov_tot5perc %>% filter(cov>= 10) %>% nrow()/ nrow(UTRCov_tot5perc)
utrCov_tot5perc_100= UTRCov_tot5perc %>% filter(cov>= 100) %>% nrow()/ nrow(UTRCov_tot5perc)
utrCov_tot5perc_1000= UTRCov_tot5perc %>% filter(cov>= 1000) %>% nrow()/ nrow(UTRCov_tot5perc)


utrCov_nuc5perc_0= UTRCov_nuc5perc %>% filter(cov> 0) %>% nrow()/ nrow(UTRCov_nuc5perc)
utrCov_nuc5perc_10= UTRCov_nuc5perc %>% filter(cov>= 10) %>% nrow()/ nrow(UTRCov_nuc5perc)
utrCov_nuc5perc_100= UTRCov_nuc5perc %>% filter(cov>= 100) %>% nrow()/ nrow(UTRCov_nuc5perc)
utrCov_nuc5perc_1000= UTRCov_nuc5perc %>% filter(cov>= 1000) %>% nrow()/ nrow(UTRCov_nuc5perc)


Total_UTR=c(utrCov_tot5perc_0, utrCov_tot5perc_10, utrCov_tot5perc_100, utrCov_tot5perc_1000)
Nuclear_UTR=c(utrCov_nuc5perc_0,utrCov_nuc5perc_10,utrCov_nuc5perc_100,utrCov_nuc5perc_1000)

covUTR_5perDF=as.data.frame(cbind(Reads,Total=Total_UTR,Nuclear=Nuclear_UTR))
covUTR_5perDF$Total=as.numeric(as.character(covUTR_5perDF$Total))
covUTR_5perDF$Nuclear=as.numeric(as.character(covUTR_5perDF$Nuclear))

covUTR_5perDF_melt=melt(covUTR_5perDF, id.vars = "Reads")
colnames(covUTR_5perDF_melt)=c("ReadCutoff", "Fraction", "Proportion" )

ggplot(covUTR_5perDF_melt, aes(x=ReadCutoff, fill=Fraction, y=Proportion, by=Fraction)) + geom_bar(position="dodge",stat="identity" ) + labs(y="Proportion of PAS covered", x="At least X reads ending at PAS", title="Proportion of PAS with PacBio read ending in it \n UTR PAS with 5% mean Usage") + scale_fill_brewer(palette = "Set1") + annotate("text", label="Total PAS = 5,984", x="1000 reads", y=.7) + annotate("text", label="Nuclear PAS = 6,793", x="1000 reads", y=.6)

Past versions of unnamed-chunk-16-1.png

Make the first plot (utr v intron for the 5%) Process in excel

allProp=read.table("../data/pacbio_cov/PacBioPropCov.txt", head=T, stringsAsFactors = F)
allProp$Read=as.factor(allProp$Read)
allPropPlot=ggplot(allProp, aes(x=Read, y=Proportion, by=Location, fill=Location)) + geom_bar(position="dodge",stat="identity" ) +facet_grid(~Fraction) + labs(y="Proportion of PAS covered", x="At least X reads ending at PAS", title="Proportion of PAS with PacBio read ending in it \n PAS with 5% mean Usage")+scale_fill_manual(values=c("red","blue"))

allPropPlot

Past versions of unnamed-chunk-17-1.png

ggsave(allPropPlot, filename = "../output/plots/PacBioReadsEndingAtPAS_5percCov.png")

Saving 7 x 5 in image

UTRCov_tot5perc %>% semi_join(UTRCov_nuc5perc, by="peak") %>% nrow()

[1] 5538

intronCov_tot5perc %>% semi_join(intronCov_nuc5perc, by="peak") %>% nrow()

[1] 15270

tot5Perc_peaks =tot5Perc %>% select(peak)
nuc5Perc_peaks=nuc5Perc %>% select(peak)
nrow(tot5Perc)

[1] 33002

nrow(nuc5Perc)

[1] 37370

ineither=tot5Perc_peaks %>% full_join(nuc5Perc_peaks, by="peak")
nrow(ineither)

[1] 40066

Pick lines for ours

I need to look at the top nuclear specific QTLs and chose individuals with both homozygous alleles.

totQTL=read.table("../data/ApaQTLs/TotalapaQTLs.GeneLocAnno.noMP.5perc.10FDR.txt", head=F,stringsAsFactors = F, col.names = c('pid', 'nvar', 'shape1', 'shape2' ,'dummy' ,'sid' ,'dist', 'npval', 'slope', 'ppval', 'bpval' ,'bh')) %>% separate(pid, into=c("chr", "start", "end", "geneID"), sep=":") %>% separate(geneID, into=c("gene", "strand", "peak"),sep="_") 

NucQTL=read.table("../data/ApaQTLs/NuclearapaQTLs.GeneLocAnno.noMP.5perc.10FDR.txt", head=F,stringsAsFactors = F, col.names = c('pid', 'nvar', 'shape1', 'shape2' ,'dummy' ,'sid' ,'dist', 'npval', 'slope', 'ppval', 'bpval' ,'bh')) %>% separate(pid, into=c("chr", "start", "end", "geneID"), sep=":") %>% separate(geneID, into=c("gene", "strand", "peak"),sep="_") 

Filter out the genes with a total qtl from the nuclear qlts genes:

NucOnlyQTL= NucQTL %>% anti_join(totQTL, by="gene") %>% arrange(slope) 

head(NucOnlyQTL)

  chr     start       end      gene strand       peak nvar   shape1
1   1 149585233 149585320 LINC00869      -   peak9493   40 1.066010
2   1 249211696 249211944     PGBD2      -  peak16578  142 0.991565
3   6  32519163  32519261  HLA-DRB6      + peak138621 2701 1.003040
4   1 113667733 113667823     LRIG2      -   peak8001  123 0.947699
5   6 109765267 109765350    MICAL1      + peak142627  122 0.886992
6  19   6535487   6535939    TNFSF9      -  peak77427  219 0.986695
    shape2   dummy         sid   dist       npval     slope       ppval
1  6.36770 48.4641 1:149601043  15809 3.31180e-04 -13.63700 0.000999001
2 27.35320 40.7121 1:249232758  21061 4.94006e-08  -7.46381 0.000999001
3  6.50355 35.2570  6:32505477 -13687 1.34996e-06  -2.94929 0.000999001
4 15.94130 42.7219 1:113678884  11150 3.66278e-09  -2.40161 0.000999001
5 15.36680 46.6074 6:109754433 -10835 1.03391e-05  -2.10134 0.001998000
6 34.27590 39.0277  19:6551945  16457 2.70575e-07  -2.05682 0.000999001
        bpval           bh
1 1.48525e-03 0.0949692158
2 2.05616e-05 0.0036973232
3 2.71955e-04 0.0297604242
4 1.27763e-06 0.0003351679
5 6.94153e-04 0.0569309509
6 1.76541e-04 0.0210644761

See if any of these are also nuclear enriched from diff iso:

nucMoreUsed=read.table("../data/diff_iso_GeneLocAnno/SigPeaksHigherInNuc.txt", header = T) %>% separate(intron, into=c("chr", "start", "end", "gene"),sep=":")

Join the nuc only by these peaks:

NucOnlyQTL_diff=NucOnlyQTL %>% inner_join(nucMoreUsed, by=c("start", "end", "gene"))%>% arrange(abs(slope))   %>% select(gene, peak, sid, slope,deltapsi )

NucOnlyQTL_diff

           gene       peak         sid     slope   deltapsi
1        IFNGR2 peak101664 21:34785672  0.460097 -0.2277754
2         RCSD1  peak11626 1:167610743  0.524358 -0.2000519
3          STX8  peak66383  17:9481747 -0.542269 -0.2115276
4         ZNRD1 peak138165  6:30017632 -0.774940 -0.2165452
5       FAM207A peak102641 21:46396659  0.923893 -0.2533298
6         DOCK7   peak5329  1:63018852  0.935357 -0.2357179
7        SMURF2  peak71606 17:62663377  1.102180 -0.3116123
8     LINC00476 peak167478  9:98613783 -1.105940 -0.2108329
9         ACTR2  peak86635  2:65478371 -1.112100 -0.2392005
10        FOXN2  peak85407  2:48547443 -1.192880 -0.2458012
11         DPF2  peak28278 11:65109233 -1.198550 -0.2090628
12        CDC42   peak1782  1:22435723  1.279150 -0.2717837
13       ZNF738  peak79187 19:21543803  1.489360 -0.2156677
14 LOC101927151  peak79354 19:28291331 -1.685500 -0.2078720
15       RPUSD2  peak53419 15:40861040 -1.716780 -0.2243680
16         NFIC  peak77091  19:3383733 -1.885840 -0.2067161

Look at these and the snps:

RPUSD2: 15:40861040
ACTR2 2:65478371 SMURF2 17:62663377 ZNF738 19:21543803

Pull in the genotypes for these:

headerGeno=read.table("../data/pacbio_cov/indiv.geno.header.txt", stringsAsFactors = F) %>% t() 

headerGeno= as.vector(headerGeno)

RPUSD2_gene=read.table("../data/pacbio_cov/geno_15:408610.txt", stringsAsFactors = F, col.names  =headerGeno)

ACTR2_gene=read.table("../data/pacbio_cov/geno_2:65478371.txt", stringsAsFactors = F, col.names  =headerGeno)

ZNF738_geno=read.table("../data/pacbio_cov/geno_19:21543803.txt", stringsAsFactors = F, col.names  =headerGeno )

SMURF2_geno=read.table("../data/pacbio_cov/geno_17:62663377.txt", stringsAsFactors = F, col.names  =headerGeno )

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] 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        pillar_1.3.0       glue_1.3.0        
[10] withr_2.1.2        RColorBrewer_1.1-2 modelr_0.1.2      
[13] readxl_1.1.0       bindr_0.1.1        plyr_1.8.4        
[16] munsell_0.5.0      gtable_0.2.0       cellranger_1.1.0  
[19] rvest_0.3.2        evaluate_0.13      labeling_0.3      
[22] knitr_1.20         broom_0.5.0        Rcpp_0.12.19      
[25] scales_1.0.0       backports_1.1.2    jsonlite_1.6      
[28] fs_1.2.6           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.11     httr_1.3.1        
[46] rstudioapi_0.9.0   R6_2.3.0           nlme_3.1-137      
[49] git2r_0.24.0       compiler_3.5.1