Compare new and old peaks

Total

totalPeakUs=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Total.fc", header = T, stringsAsFactors = F) %>% separate(chrom, sep = ":", into = c("chr", "start", "end", "id")) %>% separate(id, sep="_", into=c("gene", "loc", "strand", "peak"))

Warning: Expected 4 pieces. Additional pieces discarded in 4 rows [14735,
14736, 14737, 14738].

ind=colnames(totalPeakUs)[8:dim(totalPeakUs)[2]]
totalPeakUs_CountNum=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Total.CountsOnlyNumeric", col.names = ind)


#numeric with anno
totalPeak=as.data.frame(cbind(totalPeakUs[,1:7], totalPeakUs_CountNum))

totalPeakUs_CountNum_mean=rowMeans(totalPeakUs_CountNum)

#append mean to anno
TotalPeakUSMean=as.data.frame(cbind(totalPeakUs[,1:7],mean=totalPeakUs_CountNum_mean))

uniqnewPasnum=uniqnew  %>% separate(ID ,into=c("peaknum", "geneloc"),sep=":") %>% mutate(peak=paste("peak", peaknum, sep="")) %>% select(peak)

Filter these inthe mean usage:

TotalPeakUSMeanClass= TotalPeakUSMean %>% mutate(New=ifelse(peak %in% uniqnewPasnum$peak,"new", "original")) %>% mutate(Cutoff=ifelse(mean>=.05, "Yes", "No"))

mean(TotalPeakUSMean$mean)

[1] 0.2378282

Plot:

ggplot(TotalPeakUSMeanClass, aes(y=mean,x=New)) + geom_violin() + geom_hline(yintercept = mean(TotalPeakUSMean$mean), col="red") +  geom_hline(yintercept = .05, col="Blue")

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

TotalPeakUSMeanClass_newonly= TotalPeakUSMeanClass %>% filter(New=="new")

ggplot(TotalPeakUSMeanClass_newonly, aes(y=mean, x="")) + geom_violin() + geom_hline(yintercept = mean(TotalPeakUSMean$mean), col="red") +  geom_hline(yintercept = .05, col="Blue") + labs(x="", y="Mean Usage", title="Mean Usage of New PAS")

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

ggplot(TotalPeakUSMeanClass_newonly, aes(x=mean)) + geom_histogram(bins=50) + geom_vline(xintercept = mean(TotalPeakUSMean$mean), col="red")

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

This shows me the new peaks are the peaks that barely passed the cuttoff before.

write out file with information about new and old peaks

Peak_newOld=TotalPeakUSMeanClass %>% select(-mean)
write.table(Peak_newOld, file="../data/peaks_5perc/NewVOldPeaks.txt", col.names = T, row.names = F, quote=F)

Total:

First I will upload the usage dataframe.

totCounts=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Total.fc", stringsAsFactors = F, header = T)
ind=colnames(totCounts)[2:55]
totUsage=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Total.CountsOnlyNumeric", stringsAsFactors = F, header = F,col.names = ind)



batch1.2.3=read.table("../data/MetaDataSequencing.txt", header=T,stringsAsFactors = F)%>% filter(fraction=="total") %>%  select(line, batch) %>% filter(batch != 4)
oldind=batch1.2.3$line

batch4=read.table("../data/MetaDataSequencing.txt", header=T,stringsAsFactors = F)%>% filter(fraction=="total") %>%  select(line, batch) %>% filter(batch == 4)

newInd=batch4$line

UsageHist= function(indVec,title,totUsage=totUsage){
  totUsage_ind=totUsage %>% select(indVec)
  meanVec=rowMeans(totUsage_ind)
  hist(meanVec, main=title,xlab="Mean Usage")
}

RUn this for different itterations of individuals:

Pick 15 random individuals from old:

sampl1=sample(oldind, 15)
sampl2=sample(oldind, 15)
sampl3=sample(oldind, 15)
sampl4=sample(oldind, 15)

png("../output/newtot.png")
par(mfrow=c(3,2))
UsageHist(indVec=newInd,title="Total Usage New (15ind)",totUsage=totUsage)
UsageHist(indVec=oldind,title="Total Usage Old (39 ind)",totUsage=totUsage)
UsageHist(indVec=sampl1,title="Total Usage Sample 15 Old",totUsage=totUsage)
UsageHist(indVec=sampl2,title="Total Usage Sample 15 Old",totUsage=totUsage)
UsageHist(indVec=sampl3,title="Total Usage Sample 15 Old",totUsage=totUsage)
UsageHist(indVec=sampl4,title="Total Usage Sample 15 Old",totUsage=totUsage)
dev.off()

png 
  2 

Nuclear

nucCounts=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Nuclear.fc", stringsAsFactors = F, header = T)
ind=colnames(nucCounts)[2:55]
nucUsage=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Nuclear.CountsOnlyNumeric", stringsAsFactors = F, header = F,col.names = ind)

UsageHist_nuc= function(indVec,title,nucUsage=nucUsage){
  nucUsage_ind=nucUsage %>% select(indVec)
  meanVec=rowMeans(nucUsage_ind)
  hist(meanVec, main=title,xlab="Mean Usage")
}

png("../output/newnuc.png")
par(mfrow=c(3,2))
UsageHist_nuc(indVec=newInd,title="Nuclear Usage New (15ind)",nucUsage=nucUsage)
UsageHist_nuc(indVec=oldind,title="Nuclear Usage Old (39ind)",nucUsage=nucUsage)
UsageHist_nuc(indVec=sampl1,title="Nuclear Usage Sample 15 Old",nucUsage=nucUsage)
UsageHist_nuc(indVec=sampl2,title="Nuclear Usage Sample 15 Old",nucUsage=nucUsage)
UsageHist_nuc(indVec=sampl3,title="Nuclear Usage Sample 15 Old",nucUsage=nucUsage)
UsageHist_nuc(indVec=sampl4,title="Nuclear Usage Sample 15 Old",nucUsage=nucUsage)
dev.off()

png 
  2 

Old

oldtotalCount=read.table("../../threeprimeseq/data/phenotypes_filtPeakTranscript_noMP_GeneLocAnno/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Total.fixed.pheno.fc",header=T,stringsAsFactors = F)
indOld=colnames(oldtotalCount)[2:56]
oldtotalUsage=read.table("../../threeprimeseq/data/phenotypes_filtPeakTranscript_noMP_GeneLocAnno/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Total.fixed.pheno.CountsOnlyNumeric.txt",col.names = indOld,stringsAsFactors = F)

png("../output/oldtot.png")
par(mfrow=c(3,2))
UsageHist(indVec=newInd,title="Old total Usage (15ind)",totUsage=oldtotalUsage)
UsageHist(indVec=oldind,title="Old total Usage (39ind)",totUsage=oldtotalUsage)
UsageHist(indVec=sampl1,title="Old total Usage sample 15 ind",totUsage=oldtotalUsage)
UsageHist(indVec=sampl2,title="Old total Usage sample 15 ind",totUsage=oldtotalUsage)
UsageHist(indVec=sampl3,title="Old total Usage sample 15 ind",totUsage=oldtotalUsage)
UsageHist(indVec=sampl4,title="Old total Usage sample 15 ind",totUsage=oldtotalUsage)
dev.off()

png 
  2 

oldnuclearCount=read.table("../../threeprimeseq/data/phenotypes_filtPeakTranscript_noMP_GeneLocAnno/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Nuclear.fixed.pheno.fc",header=T,stringsAsFactors = F)
indOldN=colnames(oldnuclearCount)[2:56]
oldnuclearUsage=read.table("../../threeprimeseq/data/phenotypes_filtPeakTranscript_noMP_GeneLocAnno/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Nuclear.fixed.pheno.CountsOnlyNumeric.txt",col.names = indOldN,stringsAsFactors = F)

png("../output/oldnuc.png")
par(mfrow=c(3,2))
UsageHist_nuc(indVec=newInd,title="Old nuclear Usage (15ind)",nucUsage=oldnuclearUsage)
UsageHist_nuc(indVec=oldind,title="Old nuclear Usage (39ind)",nucUsage=oldnuclearUsage)
UsageHist_nuc(indVec=sampl1,title="Old nuclear Usage sample 15 ind",nucUsage=oldnuclearUsage)
UsageHist_nuc(indVec=sampl2,title="Old nuclear Usage sample 15 ind",nucUsage=oldnuclearUsage)
UsageHist_nuc(indVec=sampl3,title="Old nuclear Usage sample 15 ind",nucUsage=oldnuclearUsage)
UsageHist_nuc(indVec=sampl4,title="Old nuclear Usage sample 15 ind",nucUsage=oldnuclearUsage)
dev.off()

png 
  2 

Total

totalPeak_inNew=totalPeak %>% filter(peak %in% uniqnewPasnum$peak)

totalPeak_inNew_melt=melt(totalPeak_inNew,id.vars=c("chr", "start","end", "gene","loc", "peak","strand"),value.name = "Usage" ,variable.name = "Ind") %>% mutate(New15=ifelse(Ind %in% batch4$line, "Yes", "No"))


ggplot(totalPeak_inNew_melt, aes(x=New15, y=Usage, fill=New15))+geom_boxplot(width=.5) + theme(legend.position = "none") + labs(x="Individuals in New batch")

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

Look if means are different:

totalPeak_inNew_meltGroup=totalPeak_inNew_melt %>% group_by(New15,peak) %>% summarise(meanUsage=mean(Usage))

ggplot(totalPeak_inNew_meltGroup, aes(x=New15, y=meanUsage, fill=New15))+geom_boxplot(width=.5) + theme(legend.position = "none") + labs(x="Individuals in New batch", y="Mean Usage in Group", title="Mean usage for new PAS in 39 ind v new 15")

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

ggplot(totalPeak_inNew_meltGroup, aes(x=meanUsage, group=New15, fill=New15))+geom_density(alpha=.3)  + labs( x="Mean Usage in Group", title="Mean usage for new PAS in 39 ind v new 15")+ scale_fill_discrete(name = "Ind in new 15")

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

totalPeak_inNew_meltGroup_spre=totalPeak_inNew_meltGroup %>% spread(New15,meanUsage)


ggplot(totalPeak_inNew_meltGroup_spre,aes(x=No, y=Yes)) + geom_point() + geom_density_2d() + geom_smooth(method = "lm")+annotate("text", label="New=old+0.053680" ,x=.8,y=0) + labs(x="Original 39 Mean", y="Original 15 Mean", title="Mean Usage in New data by Ind set")

summary(lm(data=totalPeak_inNew_meltGroup_spre, Yes~No))

Call:
lm(formula = Yes ~ No, data = totalPeak_inNew_meltGroup_spre)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.80920 -0.08513 -0.03408  0.05478  0.86924 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.053680   0.001055   50.86   <2e-16 ***
No          1.002046   0.003416  293.35   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1584 on 42931 degrees of freedom
Multiple R-squared:  0.6672,    Adjusted R-squared:  0.6671 
F-statistic: 8.605e+04 on 1 and 42931 DF,  p-value: < 2.2e-16

Nuclear

nuclearPeakUs=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Nuclear.fc", header = T, stringsAsFactors = F) %>% separate(chrom, sep = ":", into = c("chr", "start", "end", "id")) %>% separate(id, sep="_", into=c("gene", "loc", "strand", "peak"))

Warning: Expected 4 pieces. Additional pieces discarded in 4 rows [14735,
14736, 14737, 14738].

ind=colnames(nuclearPeakUs)[8:dim(nuclearPeakUs)[2]]
nuclearPeakUs_CountNum=read.table("../data/phenotype/APApeak_Phenotype_GeneLocAnno.Nuclear.CountsOnlyNumeric", col.names = ind)


#numeric with anno
nuclearPeak=as.data.frame(cbind(nuclearPeakUs[,1:7], nuclearPeakUs_CountNum))



nuclearPeak_inNew=nuclearPeak %>% filter(peak %in% uniqnewPasnum$peak)

nuclearPeak_inNew_melt=melt(nuclearPeak_inNew,id.vars=c("chr", "start","end", "gene","loc", "peak","strand"),value.name = "Usage" ,variable.name = "Ind") %>% mutate(New15=ifelse(Ind %in% batch4$line, "Yes", "No"))

nuclearPeak_inNew_meltGroup=nuclearPeak_inNew_melt %>% group_by(New15,peak) %>% summarise(meanUsage=mean(Usage))

ggplot(nuclearPeak_inNew_meltGroup, aes(x=New15, y=meanUsage, fill=New15))+geom_boxplot(width=.5) + theme(legend.position = "none") + labs(x="Individuals in New batch", y="Mean Usage in Group", title="Nuclear Mean usage for new PAS in 39 ind v new 15")

ggplot(nuclearPeak_inNew_meltGroup, aes(x=meanUsage, group=New15, fill=New15))+geom_density(alpha=.3)  + labs( x="Mean Usage in Group", title="Nuclear Mean usage for new PAS in 39 ind v new 15")+ scale_fill_discrete(name = "Ind in new 15")