Last updated: 2020-03-18

Checks: 7 0

Knit directory: Comparative_APA/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.6.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20190902) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    code/chimp_log/
    Ignored:    code/human_log/
    Ignored:    data/.DS_Store
    Ignored:    data/TrialFiltersMeta.txt.sb-9845453e-R58Y0Q/
    Ignored:    data/mediation_prot/
    Ignored:    data/metadata_HCpanel.txt.sb-a5794dd2-i594qs/
    Ignored:    output/.DS_Store

Untracked files:
    Untracked:  ._.DS_Store
    Untracked:  Chimp/
    Untracked:  Human/
    Untracked:  analysis/CrossChimpThreePrime.Rmd
    Untracked:  analysis/DiffTransProtvsExpression.Rmd
    Untracked:  analysis/DiffUsedUTR.Rmd
    Untracked:  analysis/GvizPlots.Rmd
    Untracked:  analysis/HandC.TvN
    Untracked:  analysis/PhenotypeOverlap10.Rmd
    Untracked:  analysis/annotationBias.Rmd
    Untracked:  analysis/assessReadQual.Rmd
    Untracked:  analysis/diffExpressionPantro6.Rmd
    Untracked:  analysis/orthoexonAnno.Rmd
    Untracked:  analysis/pol2.Rmd
    Untracked:  code/._ClassifyLeafviz.sh
    Untracked:  code/._Config_chimp.yaml
    Untracked:  code/._Config_chimp_full.yaml
    Untracked:  code/._Config_human.yaml
    Untracked:  code/._ConvertJunc2Bed.sh
    Untracked:  code/._CountNucleotides.py
    Untracked:  code/._CrossMapChimpRNA.sh
    Untracked:  code/._CrossMapThreeprime.sh
    Untracked:  code/._DiffSplice.sh
    Untracked:  code/._DiffSplicePlots.sh
    Untracked:  code/._DiffSplicePlots_gencode.sh
    Untracked:  code/._DiffSplice_gencode.sh
    Untracked:  code/._DiffSplice_removebad.sh
    Untracked:  code/._FindIntronForDomPAS.sh
    Untracked:  code/._FindIntronForDomPAS_DF.sh
    Untracked:  code/._GetMAPQscore.py
    Untracked:  code/._GetSecondaryMap.py
    Untracked:  code/._Lift5perPAS.sh
    Untracked:  code/._LiftFinalChimpJunc2Human.sh
    Untracked:  code/._LiftOrthoPAS2chimp.sh
    Untracked:  code/._MapBadSamples.sh
    Untracked:  code/._PAS_ATTAAA.sh
    Untracked:  code/._PAS_ATTAAA_df.sh
    Untracked:  code/._PAS_seqExpanded.sh
    Untracked:  code/._PASsequences.sh
    Untracked:  code/._PASsequences_DF.sh
    Untracked:  code/._PlotNuclearUsagebySpecies.R
    Untracked:  code/._PlotNuclearUsagebySpecies_DF.R
    Untracked:  code/._QuantMergedClusters.sh
    Untracked:  code/._RNATranscriptDTplot.sh
    Untracked:  code/._ReverseLiftFilter.R
    Untracked:  code/._RunFixLeafCluster.sh
    Untracked:  code/._RunNegMCMediation.sh
    Untracked:  code/._RunNegMCMediationDF.sh
    Untracked:  code/._RunPosMCMediationDF.err
    Untracked:  code/._RunPosMCMediationDF.sh
    Untracked:  code/._SAF2Bed.py
    Untracked:  code/._Snakefile
    Untracked:  code/._SnakefilePAS
    Untracked:  code/._SnakefilePASfilt
    Untracked:  code/._SortIndexBadSamples.sh
    Untracked:  code/._assignPeak2Intronicregion
    Untracked:  code/._assignPeak2Intronicregion.sh
    Untracked:  code/._bed215upbed.py
    Untracked:  code/._bed2SAF_gen.py
    Untracked:  code/._buildIndecpantro5
    Untracked:  code/._buildIndecpantro5.sh
    Untracked:  code/._buildLeafviz.sh
    Untracked:  code/._buildLeafviz_leadAnno.sh
    Untracked:  code/._buildStarIndex.sh
    Untracked:  code/._chimpChromprder.sh
    Untracked:  code/._chooseSignalSite.py
    Untracked:  code/._cleanbed2saf.py
    Untracked:  code/._cluster.json
    Untracked:  code/._cluster2bed.py
    Untracked:  code/._clusterLiftReverse.sh
    Untracked:  code/._clusterLiftReverse_removebad.sh
    Untracked:  code/._clusterLiftprimary.sh
    Untracked:  code/._clusterLiftprimary_removebad.sh
    Untracked:  code/._converBam2Junc.sh
    Untracked:  code/._converBam2Junc_removeBad.sh
    Untracked:  code/._extraSnakefiltpas
    Untracked:  code/._extractPhyloReg.py
    Untracked:  code/._extractPhyloRegGene.py
    Untracked:  code/._extractPhylopGeneral.ph
    Untracked:  code/._extractPhylopGeneral.py
    Untracked:  code/._extractPhylopReg200down.py
    Untracked:  code/._extractPhylopReg200up.py
    Untracked:  code/._filter5percPAS.py
    Untracked:  code/._filterNumChroms.py
    Untracked:  code/._filterPASforMP.py
    Untracked:  code/._filterPostLift.py
    Untracked:  code/._fixExonFC.py
    Untracked:  code/._fixLeafCluster.py
    Untracked:  code/._fixLiftedJunc.py
    Untracked:  code/._fixUTRexonanno.py
    Untracked:  code/._formathg38Anno.py
    Untracked:  code/._formatpantro6Anno.py
    Untracked:  code/._getRNAseqMapStats.sh
    Untracked:  code/._hg19MapStats.sh
    Untracked:  code/._humanChromorder.sh
    Untracked:  code/._intersectLiftedPAS.sh
    Untracked:  code/._liftJunctionFiles.sh
    Untracked:  code/._liftPAS19to38.sh
    Untracked:  code/._liftedchimpJunc2human.sh
    Untracked:  code/._makeNuclearDapaplots.sh
    Untracked:  code/._makeNuclearDapaplots_DF.sh
    Untracked:  code/._makeSamplyGroupsHuman_TvN.py
    Untracked:  code/._mapRNAseqhg19.sh
    Untracked:  code/._mapRNAseqhg19_newPipeline.sh
    Untracked:  code/._maphg19.sh
    Untracked:  code/._maphg19_subjunc.sh
    Untracked:  code/._mediation_test.R
    Untracked:  code/._mergeChimp3prime_inhg38.sh
    Untracked:  code/._mergeandBWRNAseq.sh
    Untracked:  code/._mergedBam2BW.sh
    Untracked:  code/._nameClusters.py
    Untracked:  code/._negativeMediation_montecarlo.R
    Untracked:  code/._negativeMediation_montecarloDF.R
    Untracked:  code/._numMultimap.py
    Untracked:  code/._overlapapaQTLPAS.sh
    Untracked:  code/._parseHg38.py
    Untracked:  code/._postiveMediation_montecarlo_DF.R
    Untracked:  code/._prepareCleanLiftedFC_5perc4LC.py
    Untracked:  code/._prepareLeafvizAnno.sh
    Untracked:  code/._preparePAS4lift.py
    Untracked:  code/._primaryLift.sh
    Untracked:  code/._processhg38exons.py
    Untracked:  code/._quantJunc.sh
    Untracked:  code/._quantJunc_TEST.sh
    Untracked:  code/._quantJunc_removeBad.sh
    Untracked:  code/._quantMerged_seperatly.sh
    Untracked:  code/._recLiftchim2human.sh
    Untracked:  code/._revLiftPAShg38to19.sh
    Untracked:  code/._reverseLift.sh
    Untracked:  code/._runCheckReverseLift.sh
    Untracked:  code/._runChimpDiffIso.sh
    Untracked:  code/._runCountNucleotides.sh
    Untracked:  code/._runFilterNumChroms.sh
    Untracked:  code/._runHumanDiffIso.sh
    Untracked:  code/._runNuclearDiffIso_DF.sh
    Untracked:  code/._runNuclearDifffIso.sh
    Untracked:  code/._runTotalDiffIso.sh
    Untracked:  code/._run_chimpverifybam.sh
    Untracked:  code/._run_verifyBam.sh
    Untracked:  code/._snakemake.batch
    Untracked:  code/._snakemakePAS.batch
    Untracked:  code/._snakemakePASchimp.batch
    Untracked:  code/._snakemakePAShuman.batch
    Untracked:  code/._snakemake_chimp.batch
    Untracked:  code/._snakemake_human.batch
    Untracked:  code/._snakemakefiltPAS.batch
    Untracked:  code/._snakemakefiltPAS_chimp
    Untracked:  code/._snakemakefiltPAS_chimp.sh
    Untracked:  code/._snakemakefiltPAS_human.sh
    Untracked:  code/._spliceSite2Fasta.py
    Untracked:  code/._submit-snakemake-chimp.sh
    Untracked:  code/._submit-snakemake-human.sh
    Untracked:  code/._submit-snakemakePAS-chimp.sh
    Untracked:  code/._submit-snakemakePAS-human.sh
    Untracked:  code/._submit-snakemakefiltPAS-chimp.sh
    Untracked:  code/._submit-snakemakefiltPAS-human.sh
    Untracked:  code/._subset_diffisopheno_Nuclear_HvC.py
    Untracked:  code/._subset_diffisopheno_Nuclear_HvC_DF.py
    Untracked:  code/._subset_diffisopheno_Total_HvC.py
    Untracked:  code/._threeprimeOrthoFC.sh
    Untracked:  code/._transcriptDTplotsNuclear.sh
    Untracked:  code/._verifyBam4973.sh
    Untracked:  code/._verifyBam4973inHuman.sh
    Untracked:  code/._wrap_chimpverifybam.sh
    Untracked:  code/._wrap_verifyBam.sh
    Untracked:  code/._writeMergecode.py
    Untracked:  code/.snakemake/
    Untracked:  code/ClassifyLeafviz.sh
    Untracked:  code/Config_chimp.yaml
    Untracked:  code/Config_chimp_full.yaml
    Untracked:  code/Config_human.yaml
    Untracked:  code/ConvertJunc2Bed.err
    Untracked:  code/ConvertJunc2Bed.out
    Untracked:  code/ConvertJunc2Bed.sh
    Untracked:  code/CountNucleotides.py
    Untracked:  code/CrossMapChimpRNA.sh
    Untracked:  code/CrossMapThreeprime.sh
    Untracked:  code/CrossmapChimp3prime.err
    Untracked:  code/CrossmapChimp3prime.out
    Untracked:  code/CrossmapChimpRNA.err
    Untracked:  code/CrossmapChimpRNA.out
    Untracked:  code/DiffSplice.err
    Untracked:  code/DiffSplice.out
    Untracked:  code/DiffSplice.sh
    Untracked:  code/DiffSplicePlots.err
    Untracked:  code/DiffSplicePlots.out
    Untracked:  code/DiffSplicePlots.sh
    Untracked:  code/DiffSplicePlots_gencode.sh
    Untracked:  code/DiffSplice_gencode.sh
    Untracked:  code/DiffSplice_removebad.err
    Untracked:  code/DiffSplice_removebad.out
    Untracked:  code/DiffSplice_removebad.sh
    Untracked:  code/FilterReverseLift.err
    Untracked:  code/FilterReverseLift.out
    Untracked:  code/FindIntronForDomPAS.err
    Untracked:  code/FindIntronForDomPAS.out
    Untracked:  code/FindIntronForDomPAS.sh
    Untracked:  code/FindIntronForDomPAS_DF.sh
    Untracked:  code/GencodeDiffSplice.err
    Untracked:  code/GencodeDiffSplice.out
    Untracked:  code/GetMAPQscore.py
    Untracked:  code/GetSecondaryMap.py
    Untracked:  code/HchromOrder.err
    Untracked:  code/HchromOrder.out
    Untracked:  code/JunctionLift.err
    Untracked:  code/JunctionLift.out
    Untracked:  code/JunctionLiftFinalChimp.err
    Untracked:  code/JunctionLiftFinalChimp.out
    Untracked:  code/Lift5perPAS.sh
    Untracked:  code/Lift5perPASbed.err
    Untracked:  code/Lift5perPASbed.out
    Untracked:  code/LiftClustersFirst.err
    Untracked:  code/LiftClustersFirst.out
    Untracked:  code/LiftClustersFirst_remove.err
    Untracked:  code/LiftClustersFirst_remove.out
    Untracked:  code/LiftClustersSecond.err
    Untracked:  code/LiftClustersSecond.out
    Untracked:  code/LiftClustersSecond_remove.err
    Untracked:  code/LiftClustersSecond_remove.out
    Untracked:  code/LiftFinalChimpJunc2Human.sh
    Untracked:  code/LiftOrthoPAS2chimp.sh
    Untracked:  code/LiftorthoPAS.err
    Untracked:  code/LiftorthoPASt.out
    Untracked:  code/Log.out
    Untracked:  code/MapBadSamples.err
    Untracked:  code/MapBadSamples.out
    Untracked:  code/MapBadSamples.sh
    Untracked:  code/MapStats.err
    Untracked:  code/MapStats.out
    Untracked:  code/MaxEntCode/
    Untracked:  code/MergeClusters.err
    Untracked:  code/MergeClusters.out
    Untracked:  code/MergeClusters.sh
    Untracked:  code/PAS_ATTAAA.err
    Untracked:  code/PAS_ATTAAA.out
    Untracked:  code/PAS_ATTAAA.sh
    Untracked:  code/PAS_ATTAAADF.err
    Untracked:  code/PAS_ATTAAADF.out
    Untracked:  code/PAS_ATTAAA_df.sh
    Untracked:  code/PAS_seqExpanded.sh
    Untracked:  code/PAS_sequence.err
    Untracked:  code/PAS_sequence.out
    Untracked:  code/PAS_sequenceDF.err
    Untracked:  code/PAS_sequenceDF.out
    Untracked:  code/PASexpanded_sequenceDF.err
    Untracked:  code/PASexpanded_sequenceDF.out
    Untracked:  code/PASsequences.sh
    Untracked:  code/PASsequences_DF.sh
    Untracked:  code/PlotNuclearUsagebySpecies.R
    Untracked:  code/PlotNuclearUsagebySpecies_DF.R
    Untracked:  code/QuantMergeClusters
    Untracked:  code/QuantMergeClusters.err
    Untracked:  code/QuantMergeClusters.out
    Untracked:  code/QuantMergedClusters.sh
    Untracked:  code/RNATranscriptDTplot.err
    Untracked:  code/RNATranscriptDTplot.out
    Untracked:  code/RNATranscriptDTplot.sh
    Untracked:  code/Rev_liftoverPAShg19to38.err
    Untracked:  code/Rev_liftoverPAShg19to38.out
    Untracked:  code/ReverseLiftFilter.R
    Untracked:  code/RunFixCluster.err
    Untracked:  code/RunFixCluster.out
    Untracked:  code/RunFixLeafCluster.sh
    Untracked:  code/RunNegMCMediation.err
    Untracked:  code/RunNegMCMediation.sh
    Untracked:  code/RunNegMCMediationDF.err
    Untracked:  code/RunNegMCMediationDF.out
    Untracked:  code/RunNegMCMediationDF.sh
    Untracked:  code/RunNegMCMediationr.out
    Untracked:  code/RunPosMCMediation.err
    Untracked:  code/RunPosMCMediation.sh
    Untracked:  code/RunPosMCMediationDF.err
    Untracked:  code/RunPosMCMediationDF.out
    Untracked:  code/RunPosMCMediationDF.sh
    Untracked:  code/RunPosMCMediationr.out
    Untracked:  code/SAF215upbed_gen.py
    Untracked:  code/SAF2Bed.py
    Untracked:  code/Snakefile
    Untracked:  code/SnakefilePAS
    Untracked:  code/SnakefilePASfilt
    Untracked:  code/SortIndexBadSamples.err
    Untracked:  code/SortIndexBadSamples.out
    Untracked:  code/SortIndexBadSamples.sh
    Untracked:  code/TotalTranscriptDTplot.err
    Untracked:  code/TotalTranscriptDTplot.out
    Untracked:  code/Upstream10Bases_general.py
    Untracked:  code/apaQTLsnake.err
    Untracked:  code/apaQTLsnake.out
    Untracked:  code/apaQTLsnakePAS.err
    Untracked:  code/apaQTLsnakePAS.out
    Untracked:  code/apaQTLsnakePAShuman.err
    Untracked:  code/assignPeak2Intronicregion.err
    Untracked:  code/assignPeak2Intronicregion.out
    Untracked:  code/assignPeak2Intronicregion.sh
    Untracked:  code/bam2junc.err
    Untracked:  code/bam2junc.out
    Untracked:  code/bam2junc_remove.err
    Untracked:  code/bam2junc_remove.out
    Untracked:  code/bed215upbed.py
    Untracked:  code/bed2SAF_gen.py
    Untracked:  code/bed2saf.py
    Untracked:  code/bg_to_cov.py
    Untracked:  code/buildIndecpantro5
    Untracked:  code/buildIndecpantro5.sh
    Untracked:  code/buildLeafviz.err
    Untracked:  code/buildLeafviz.out
    Untracked:  code/buildLeafviz.sh
    Untracked:  code/buildLeafviz_leadAnno.sh
    Untracked:  code/buildLeafviz_leafanno.err
    Untracked:  code/buildLeafviz_leafanno.out
    Untracked:  code/buildStarIndex.sh
    Untracked:  code/callPeaksYL.py
    Untracked:  code/chimpChromprder.sh
    Untracked:  code/chooseAnno2Bed.py
    Untracked:  code/chooseAnno2SAF.py
    Untracked:  code/chooseSignalSite.py
    Untracked:  code/chromOrder.err
    Untracked:  code/chromOrder.out
    Untracked:  code/classifyLeafviz.err
    Untracked:  code/classifyLeafviz.out
    Untracked:  code/cleanbed2saf.py
    Untracked:  code/cluster.json
    Untracked:  code/cluster2bed.py
    Untracked:  code/clusterLiftReverse.sh
    Untracked:  code/clusterLiftReverse_removebad.sh
    Untracked:  code/clusterLiftprimary.sh
    Untracked:  code/clusterLiftprimary_removebad.sh
    Untracked:  code/clusterPAS.json
    Untracked:  code/clusterfiltPAS.json
    Untracked:  code/comands2Mege.sh
    Untracked:  code/converBam2Junc.sh
    Untracked:  code/converBam2Junc_removeBad.sh
    Untracked:  code/convertNumeric.py
    Untracked:  code/environment.yaml
    Untracked:  code/extraSnakefiltpas
    Untracked:  code/extractPhyloReg.py
    Untracked:  code/extractPhyloRegGene.py
    Untracked:  code/extractPhylopGeneral.py
    Untracked:  code/extractPhylopReg200down.py
    Untracked:  code/extractPhylopReg200up.py
    Untracked:  code/filter5perc.R
    Untracked:  code/filter5percPAS.py
    Untracked:  code/filter5percPheno.py
    Untracked:  code/filterBamforMP.pysam2_gen.py
    Untracked:  code/filterJuncChroms.err
    Untracked:  code/filterJuncChroms.out
    Untracked:  code/filterMissprimingInNuc10_gen.py
    Untracked:  code/filterNumChroms.py
    Untracked:  code/filterPASforMP.py
    Untracked:  code/filterPostLift.py
    Untracked:  code/filterSAFforMP_gen.py
    Untracked:  code/filterSortBedbyCleanedBed_gen.R
    Untracked:  code/filterpeaks.py
    Untracked:  code/fixExonFC.py
    Untracked:  code/fixFChead.py
    Untracked:  code/fixFChead_bothfrac.py
    Untracked:  code/fixLeafCluster.py
    Untracked:  code/fixLiftedJunc.py
    Untracked:  code/fixUTRexonanno.py
    Untracked:  code/formathg38Anno.py
    Untracked:  code/generateStarIndex.err
    Untracked:  code/generateStarIndex.out
    Untracked:  code/generateStarIndexHuman.err
    Untracked:  code/generateStarIndexHuman.out
    Untracked:  code/getRNAseqMapStats.sh
    Untracked:  code/hg19MapStats.err
    Untracked:  code/hg19MapStats.out
    Untracked:  code/hg19MapStats.sh
    Untracked:  code/humanChromorder.sh
    Untracked:  code/humanFiles
    Untracked:  code/intersectAnno.err
    Untracked:  code/intersectAnno.out
    Untracked:  code/intersectAnnoExt.err
    Untracked:  code/intersectAnnoExt.out
    Untracked:  code/intersectLiftedPAS.sh
    Untracked:  code/leafcutter_merge_regtools_redo.py
    Untracked:  code/liftJunctionFiles.sh
    Untracked:  code/liftPAS19to38.sh
    Untracked:  code/liftoverPAShg19to38.err
    Untracked:  code/liftoverPAShg19to38.out
    Untracked:  code/log/
    Untracked:  code/make5percPeakbed.py
    Untracked:  code/makeFileID.py
    Untracked:  code/makeNuclearDapaplots.sh
    Untracked:  code/makeNuclearDapaplots_DF.sh
    Untracked:  code/makeNuclearPlots.err
    Untracked:  code/makeNuclearPlots.out
    Untracked:  code/makeNuclearPlotsDF.err
    Untracked:  code/makeNuclearPlotsDF.out
    Untracked:  code/makePheno.py
    Untracked:  code/makeSamplyGroupsChimp_TvN.py
    Untracked:  code/makeSamplyGroupsHuman_TvN.py
    Untracked:  code/mapRNAseqhg19.sh
    Untracked:  code/mapRNAseqhg19_newPipeline.sh
    Untracked:  code/maphg19.err
    Untracked:  code/maphg19.out
    Untracked:  code/maphg19.sh
    Untracked:  code/maphg19_new.err
    Untracked:  code/maphg19_new.out
    Untracked:  code/maphg19_sub.err
    Untracked:  code/maphg19_sub.out
    Untracked:  code/maphg19_subjunc.sh
    Untracked:  code/mediation_test.R
    Untracked:  code/merge.err
    Untracked:  code/mergeChimp3prime_inhg38.sh
    Untracked:  code/merge_leafcutter_clusters_redo.py
    Untracked:  code/mergeandBWRNAseq.sh
    Untracked:  code/mergeandsort_ChimpinHuman.err
    Untracked:  code/mergeandsort_ChimpinHuman.out
    Untracked:  code/mergedBam2BW.sh
    Untracked:  code/mergedbam2bw.err
    Untracked:  code/mergedbam2bw.out
    Untracked:  code/mergedbamRNAand2bw.err
    Untracked:  code/mergedbamRNAand2bw.out
    Untracked:  code/nameClusters.py
    Untracked:  code/namePeaks.py
    Untracked:  code/negativeMediation_montecarlo.R
    Untracked:  code/negativeMediation_montecarloDF.R
    Untracked:  code/nuclearTranscriptDTplot.err
    Untracked:  code/nuclearTranscriptDTplot.out
    Untracked:  code/numMultimap.py
    Untracked:  code/overlapPAS.err
    Untracked:  code/overlapPAS.out
    Untracked:  code/overlapapaQTLPAS.sh
    Untracked:  code/overlapapaQTLPAS_extended.sh
    Untracked:  code/overlapapaQTLPAS_samples.sh
    Untracked:  code/parseHg38.py
    Untracked:  code/peak2PAS.py
    Untracked:  code/pheno2countonly.R
    Untracked:  code/postiveMediation_montecarlo.R
    Untracked:  code/postiveMediation_montecarlo_DF.R
    Untracked:  code/prepareAnnoLeafviz.err
    Untracked:  code/prepareAnnoLeafviz.out
    Untracked:  code/prepareCleanLiftedFC_5perc4LC.py
    Untracked:  code/prepareLeafvizAnno.sh
    Untracked:  code/preparePAS4lift.py
    Untracked:  code/prepare_phenotype_table.py
    Untracked:  code/primaryLift.err
    Untracked:  code/primaryLift.out
    Untracked:  code/primaryLift.sh
    Untracked:  code/processhg38exons.py
    Untracked:  code/quantJunc.sh
    Untracked:  code/quantJunc_TEST.sh
    Untracked:  code/quantJunc_removeBad.sh
    Untracked:  code/quantLiftedPAS.err
    Untracked:  code/quantLiftedPAS.out
    Untracked:  code/quantLiftedPAS.sh
    Untracked:  code/quatJunc.err
    Untracked:  code/quatJunc.out
    Untracked:  code/recChimpback2Human.err
    Untracked:  code/recChimpback2Human.out
    Untracked:  code/recLiftchim2human.sh
    Untracked:  code/revLift.err
    Untracked:  code/revLift.out
    Untracked:  code/revLiftPAShg38to19.sh
    Untracked:  code/reverseLift.sh
    Untracked:  code/runCheckReverseLift.sh
    Untracked:  code/runChimpDiffIso.sh
    Untracked:  code/runChimpDiffIsoDF.sh
    Untracked:  code/runCountNucleotides.err
    Untracked:  code/runCountNucleotides.out
    Untracked:  code/runCountNucleotides.sh
    Untracked:  code/runCountNucleotidesPantro6.err
    Untracked:  code/runCountNucleotidesPantro6.out
    Untracked:  code/runCountNucleotides_pantro6.sh
    Untracked:  code/runFilterNumChroms.sh
    Untracked:  code/runHumanDiffIso.sh
    Untracked:  code/runHumanDiffIsoDF.sh
    Untracked:  code/runNuclearDiffIso_DF.sh
    Untracked:  code/runNuclearDifffIso.sh
    Untracked:  code/runTotalDiffIso.sh
    Untracked:  code/run_Chimpleafcutter_ds.err
    Untracked:  code/run_Chimpleafcutter_ds.out
    Untracked:  code/run_Chimpverifybam.err
    Untracked:  code/run_Chimpverifybam.out
    Untracked:  code/run_Humanleafcutter_dF.err
    Untracked:  code/run_Humanleafcutter_dF.out
    Untracked:  code/run_Humanleafcutter_ds.err
    Untracked:  code/run_Humanleafcutter_ds.out
    Untracked:  code/run_Nuclearleafcutter_ds.err
    Untracked:  code/run_Nuclearleafcutter_ds.out
    Untracked:  code/run_Nuclearleafcutter_dsDF.err
    Untracked:  code/run_Nuclearleafcutter_dsDF.out
    Untracked:  code/run_Totalleafcutter_ds.err
    Untracked:  code/run_Totalleafcutter_ds.out
    Untracked:  code/run_chimpverifybam.sh
    Untracked:  code/run_verifyBam.sh
    Untracked:  code/run_verifybam.err
    Untracked:  code/run_verifybam.out
    Untracked:  code/slurm-62824013.out
    Untracked:  code/slurm-62825841.out
    Untracked:  code/slurm-62826116.out
    Untracked:  code/slurm-64108209.out
    Untracked:  code/slurm-64108521.out
    Untracked:  code/slurm-64108557.out
    Untracked:  code/snakePASChimp.err
    Untracked:  code/snakePASChimp.out
    Untracked:  code/snakePAShuman.out
    Untracked:  code/snakemake.batch
    Untracked:  code/snakemakeChimp.err
    Untracked:  code/snakemakeChimp.out
    Untracked:  code/snakemakeHuman.err
    Untracked:  code/snakemakeHuman.out
    Untracked:  code/snakemakePAS.batch
    Untracked:  code/snakemakePASFiltChimp.err
    Untracked:  code/snakemakePASFiltChimp.out
    Untracked:  code/snakemakePASFiltHuman.err
    Untracked:  code/snakemakePASFiltHuman.out
    Untracked:  code/snakemakePASchimp.batch
    Untracked:  code/snakemakePAShuman.batch
    Untracked:  code/snakemake_chimp.batch
    Untracked:  code/snakemake_human.batch
    Untracked:  code/snakemakefiltPAS.batch
    Untracked:  code/snakemakefiltPAS_chimp.sh
    Untracked:  code/snakemakefiltPAS_human.sh
    Untracked:  code/spliceSite2Fasta.py
    Untracked:  code/submit-snakemake-chimp.sh
    Untracked:  code/submit-snakemake-human.sh
    Untracked:  code/submit-snakemakePAS-chimp.sh
    Untracked:  code/submit-snakemakePAS-human.sh
    Untracked:  code/submit-snakemakefiltPAS-chimp.sh
    Untracked:  code/submit-snakemakefiltPAS-human.sh
    Untracked:  code/subset_diffisopheno.py
    Untracked:  code/subset_diffisopheno_Chimp_tvN.py
    Untracked:  code/subset_diffisopheno_Chimp_tvN_DF.py
    Untracked:  code/subset_diffisopheno_Huma_tvN.py
    Untracked:  code/subset_diffisopheno_Huma_tvN_DF.py
    Untracked:  code/subset_diffisopheno_Nuclear_HvC.py
    Untracked:  code/subset_diffisopheno_Nuclear_HvC_DF.py
    Untracked:  code/subset_diffisopheno_Total_HvC.py
    Untracked:  code/test
    Untracked:  code/test.txt
    Untracked:  code/threeprimeOrthoFC.out
    Untracked:  code/threeprimeOrthoFC.sh
    Untracked:  code/threeprimeOrthoFCcd.err
    Untracked:  code/transcriptDTplotsNuclear.sh
    Untracked:  code/transcriptDTplotsTotal.sh
    Untracked:  code/verifyBam4973.sh
    Untracked:  code/verifyBam4973inHuman.sh
    Untracked:  code/verifybam4973.err
    Untracked:  code/verifybam4973.out
    Untracked:  code/verifybam4973HumanMap.err
    Untracked:  code/verifybam4973HumanMap.out
    Untracked:  code/wrap_Chimpverifybam.err
    Untracked:  code/wrap_Chimpverifybam.out
    Untracked:  code/wrap_chimpverifybam.sh
    Untracked:  code/wrap_verifyBam.sh
    Untracked:  code/wrap_verifybam.err
    Untracked:  code/wrap_verifybam.out
    Untracked:  code/writeMergecode.py
    Untracked:  data/._.DS_Store
    Untracked:  data/._HC_filenames.txt
    Untracked:  data/._HC_filenames.txt.sb-4426323c-IKIs0S
    Untracked:  data/._HC_filenames.xlsx
    Untracked:  data/._MapPantro6_meta.txt
    Untracked:  data/._MapPantro6_meta.txt.sb-a5794dd2-Cskmlm
    Untracked:  data/._MapPantro6_meta.xlsx
    Untracked:  data/._OppositeSpeciesMap.txt
    Untracked:  data/._OppositeSpeciesMap.txt.sb-a5794dd2-mayWJf
    Untracked:  data/._OppositeSpeciesMap.xlsx
    Untracked:  data/._RNASEQ_metadata.txt
    Untracked:  data/._RNASEQ_metadata.txt.sb-4426323c-TE4ns3
    Untracked:  data/._RNASEQ_metadata.txt.sb-51f67ae1-HXp7Gq
    Untracked:  data/._RNASEQ_metadata_2Removed.txt
    Untracked:  data/._RNASEQ_metadata_2Removed.txt.sb-4426323c-a4lBwx
    Untracked:  data/._RNASEQ_metadata_2Removed.xlsx
    Untracked:  data/._RNASEQ_metadata_stranded.txt
    Untracked:  data/._RNASEQ_metadata_stranded.txt.sb-a5794dd2-D659m2
    Untracked:  data/._RNASEQ_metadata_stranded.txt.sb-a5794dd2-ImNMoY
    Untracked:  data/._RNASEQ_metadata_stranded.txt.sb-e4bf31f0-ZGnGgl
    Untracked:  data/._RNASEQ_metadata_stranded.xlsx
    Untracked:  data/._TrialFiltersMeta.txt
    Untracked:  data/._TrialFiltersMeta.txt.sb-9845453e-R58Y0Q
    Untracked:  data/._metadata_HCpanel.txt
    Untracked:  data/._metadata_HCpanel.txt.sb-a3d92a2d-b9cYoF
    Untracked:  data/._metadata_HCpanel.txt.sb-a5794dd2-i594qs
    Untracked:  data/._metadata_HCpanel.txt.sb-f4823d1e-qihGek
    Untracked:  data/._metadata_HCpanel.xlsx
    Untracked:  data/._metadata_HCpanel_frompantro5.xlsx
    Untracked:  data/._~$RNASEQ_metadata.xlsx
    Untracked:  data/._~$metadata_HCpanel.xlsx
    Untracked:  data/._.xlsx
    Untracked:  data/BaseComp/
    Untracked:  data/CompapaQTLpas/
    Untracked:  data/DNDS/
    Untracked:  data/DTmatrix/
    Untracked:  data/DiffExpression/
    Untracked:  data/DiffIso_Nuclear/
    Untracked:  data/DiffIso_Nuclear_DF/
    Untracked:  data/DiffIso_Total/
    Untracked:  data/DiffSplice/
    Untracked:  data/DiffSplice_liftedJunc/
    Untracked:  data/DiffSplice_removeBad/
    Untracked:  data/DominantPAS/
    Untracked:  data/DominantPAS_DF/
    Untracked:  data/EvalPantro5/
    Untracked:  data/HC_filenames.txt
    Untracked:  data/HC_filenames.xlsx
    Untracked:  data/Khan_prot/
    Untracked:  data/Li_eqtls/
    Untracked:  data/MapPantro6_meta.txt
    Untracked:  data/MapPantro6_meta.xlsx
    Untracked:  data/MapStats/
    Untracked:  data/NormalizedClusters/
    Untracked:  data/NuclearHvC/
    Untracked:  data/NuclearHvC_DF/
    Untracked:  data/OppositeSpeciesMap.txt
    Untracked:  data/OppositeSpeciesMap.xlsx
    Untracked:  data/OrthoExonBed/
    Untracked:  data/OverlapBenchmark/
    Untracked:  data/OverlappingPAS/
    Untracked:  data/PAS/
    Untracked:  data/PAS_doubleFilter/
    Untracked:  data/Peaks_5perc/
    Untracked:  data/Pheno_5perc/
    Untracked:  data/Pheno_5perc_DF_nuclear/
    Untracked:  data/Pheno_5perc_nuclear/
    Untracked:  data/Pheno_5perc_nuclear_old/
    Untracked:  data/Pheno_5perc_total/
    Untracked:  data/PhyloP/
    Untracked:  data/RNASEQ_metadata.txt
    Untracked:  data/RNASEQ_metadata_2Removed.txt
    Untracked:  data/RNASEQ_metadata_2Removed.xlsx
    Untracked:  data/RNASEQ_metadata_stranded.txt
    Untracked:  data/RNASEQ_metadata_stranded.txt.sb-e4bf31f0-ZGnGgl/
    Untracked:  data/RNASEQ_metadata_stranded.xlsx
    Untracked:  data/SignalSites/
    Untracked:  data/SignalSites_doublefilter/
    Untracked:  data/SpliceSite/
    Untracked:  data/Threeprime2Ortho/
    Untracked:  data/TotalHvC/
    Untracked:  data/TrialFiltersMeta.txt
    Untracked:  data/TwoBadSampleAnalysis/
    Untracked:  data/Wang_ribo/
    Untracked:  data/apaQTLGenes/
    Untracked:  data/bioGRID/
    Untracked:  data/chainFiles/
    Untracked:  data/cleanPeaks_anno/
    Untracked:  data/cleanPeaks_byspecies/
    Untracked:  data/cleanPeaks_lifted/
    Untracked:  data/files4viz_nuclear/
    Untracked:  data/files4viz_nuclear_DF/
    Untracked:  data/gviz/
    Untracked:  data/leafviz/
    Untracked:  data/liftover_files/
    Untracked:  data/mediation/
    Untracked:  data/mediation_DF/
    Untracked:  data/metadata_HCpanel.txt
    Untracked:  data/metadata_HCpanel.xlsx
    Untracked:  data/metadata_HCpanel_frompantro5.txt
    Untracked:  data/metadata_HCpanel_frompantro5.xlsx
    Untracked:  data/primaryLift/
    Untracked:  data/reverseLift/
    Untracked:  data/~$RNASEQ_metadata.xlsx
    Untracked:  data/~$metadata_HCpanel.xlsx
    Untracked:  data/.xlsx
    Untracked:  output/._.DS_Store
    Untracked:  output/dtPlots/
    Untracked:  projectNotes.Rmd
    Untracked:  proteinModelSet.Rmd

Unstaged changes:
    Modified:   analysis/DiffUsedIntronic.Rmd
    Modified:   analysis/ExploredAPA.Rmd
    Modified:   analysis/ExploredAPA_DF.Rmd
    Modified:   analysis/OppositeMap.Rmd
    Modified:   analysis/SpliceSiteStrength.Rmd
    Modified:   analysis/TotalVNuclearBothSpecies.Rmd
    Modified:   analysis/annotationInfo.Rmd
    Modified:   analysis/changeMisprimcut.Rmd
    Modified:   analysis/comp2apaQTLPAS.Rmd
    Modified:   analysis/correlationPhenos.Rmd
    Modified:   analysis/dAPAandapaQTL_DF.Rmd
    Modified:   analysis/establishCutoffs.Rmd
    Modified:   analysis/investigatePantro5.Rmd
    Modified:   analysis/multiMap.Rmd
    Modified:   analysis/speciesSpecific.Rmd
    Modified:   analysis/speciesSpecific_DF.Rmd
    Modified:   analysis/upsetter_DF.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File Version Author Date Message
Rmd d481350 brimittleman 2020-03-18 add more upand down controll
html 0bbaafe brimittleman 2020-03-05 Build site.
Rmd 0665dda brimittleman 2020-03-05 add downstream control
html 7b73ce7 brimittleman 2020-03-04 Build site.
Rmd ad95271 brimittleman 2020-03-04 add upstream control

In my initial exploration of dAPA PAS I saw they are enriched for negative phylop scores. I will explore this trend further here. I will see if intron location explain the differences.

library(tidyverse)
── Attaching packages ────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.1.1       ✔ purrr   0.3.2  
✔ tibble  2.1.1       ✔ dplyr   0.8.0.1
✔ tidyr   0.8.3       ✔ stringr 1.3.1  
✔ readr   1.3.1       ✔ forcats 0.3.0  
── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
library(ggpubr)
Loading required package: magrittr

Attaching package: 'magrittr'
The following object is masked from 'package:purrr':

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

    extract
library(reshape2)

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

    smiths
DiffUsage=read.table("../data/DiffIso_Nuclear_DF/SignifianceEitherPAS_2_Nuclear.txt", header = T, stringsAsFactors = F)

PASMeta=read.table("../data/PAS_doubleFilter/PAS_10perc_either_HumanCoord_BothUsage_meta_doubleFilter.txt", header = T, stringsAsFactors = F) %>% dplyr::select(PAS, chr, start,end, gene, loc)

DiffUsagePAS=DiffUsage %>% inner_join(PASMeta, by=c("gene","chr", "start", "end"))
phylores=read.table("../data/PhyloP/PAS_phyloP.txt", col.names = c("chr","start","end", "phyloP"), stringsAsFactors = F) %>% drop_na()
NucReswPhy=read.table("../data/DiffIso_Nuclear_DF/AllPAS_withGeneSig.txt", header = T, stringsAsFactors = F) %>% inner_join(phylores, by=c("chr","start","end"))

40756 have results of the 40776

ggplot(NucReswPhy,aes(y=phyloP, x=SigPAU2,fill=SigPAU2)) + geom_boxplot() + stat_compare_means()+ scale_fill_brewer(palette = "Dark2", name="Signficant")

Version Author Date
7b73ce7 brimittleman 2020-03-04
ggplot(NucReswPhy,aes(x=phyloP, by=SigPAU2, fill=SigPAU2)) + geom_density(alpha=.5) + scale_fill_brewer(palette = "Dark2", name="Signficant PAS") + labs(title="Mean PhyloP scores for tested PAS") + annotate("text",label="Wilcoxan, p=1.4e -5",x=6,y=.75)

Version Author Date
7b73ce7 brimittleman 2020-03-04

The significant PAS have on average lower phyloP scores.

Positive scores — Measure conservation, which is slower evolution than expected, at sites that are predicted to be conserved. Negative scores — Measure acceleration, which is faster evolution than expected, at sites that are predicted to be fast-evolving.

I can look at those with negative values:

x=nrow(NucReswPhy %>% filter(SigPAU2=="Yes", phyloP<0))
m= nrow(NucReswPhy %>% filter(phyloP<0))
n=nrow(NucReswPhy %>% filter(phyloP>=0))
k=nrow(NucReswPhy %>% filter(SigPAU2=="Yes"))


#expected
which(grepl(max(dhyper(1:x, m, n, k)), dhyper(1:x, m, n, k)))
[1] 707
#actual:
x
[1] 788
#pval
phyper(x,m,n,k,lower.tail=F)
[1] 0.0001570509
b=nrow(NucReswPhy %>% filter(SigPAU2=="Yes", phyloP<0))
n=nrow(NucReswPhy %>% filter(SigPAU2=="Yes"))
B=nrow(NucReswPhy %>% filter(phyloP<0))
N=nrow(NucReswPhy)

(b/n)/(B/N)
[1] 1.114898

This means these regions are more likely to be fast evolving.

Look at this by location: (is it driven by region)

NucReswPhy_meta= NucReswPhy %>% inner_join(PASMeta, by=c("chr", "start", "end", "gene"))

ggplot(NucReswPhy_meta,aes(x=phyloP, by=SigPAU2, fill=SigPAU2)) + geom_density(alpha=.5) + scale_fill_brewer(palette = "Dark2") + facet_grid(~loc)

Version Author Date
7b73ce7 brimittleman 2020-03-04
NucReswPhy_meta_group=NucReswPhy_meta %>% group_by(loc,SigPAU2) %>% summarise(n=n(),meanPhylo=mean(phyloP))
NucReswPhy_meta_group
# A tibble: 10 x 4
# Groups:   loc [5]
   loc    SigPAU2     n meanPhylo
   <chr>  <chr>   <int>     <dbl>
 1 cds    No       7141    2.16  
 2 cds    Yes       333    2.16  
 3 end    No       3564    0.450 
 4 end    Yes       247    0.403 
 5 intron No      10478    0.0630
 6 intron Yes       737    0.0702
 7 utr3   No      15351    1.04  
 8 utr3   Yes      1659    0.933 
 9 utr5   No       1151    0.300 
10 utr5   Yes        95    0.230 

Control sequence:upstream

(upstream 200)

Look at the 200 basepairs upstream of each PAS as a control.

metaStrand=read.table("../data/PAS_doubleFilter/PAS_10perc_either_HumanCoord_BothUsage_meta_doubleFilter.txt", header = T, stringsAsFactors = F) %>% select(chr, start,end, strandFix, PAS)

NucReswPhy_upstream=NucReswPhy %>% inner_join(metaStrand,by=c("chr", "start", "end")) %>% mutate(newStart=ifelse(strandFix=="+", start - 200, end), newEnd=ifelse(strandFix=="+", start, end +200))

NucReswPhy_upstreambed=NucReswPhy_upstream %>% select(chr, newStart, newEnd, PAS, Human, strandFix)

write.table(NucReswPhy_upstreambed,"../data/PhyloP/PAS_200upregions.bed",col.names = F,row.names = F,quote = F,sep="\t")
python extractPhylopReg200up.py
Phylo200UpContron=read.table("../data/PhyloP/PAS_phyloP_200upstream.txt",stringsAsFactors = F, col.names = c("chr", "start","end", "PAS","UpstreamControl_Phylop")) 

NucReswPhyandC=read.table("../data/DiffIso_Nuclear_DF/AllPAS_withGeneSig.txt", header = T, stringsAsFactors = F) %>% inner_join(phylores, by=c("chr", "start","end")) %>% inner_join(metaStrand,by=c("chr", "start", "end"))%>% inner_join(Phylo200UpContron, by="PAS")  %>% drop_na()

NucReswPhyandCsmall=NucReswPhyandC %>% select(PAS,SigPAU2,phyloP ,UpstreamControl_Phylop ) %>% gather("set", "Phylop", -PAS, -SigPAU2)

wilcox.test(NucReswPhyandC$phyloP, NucReswPhyandC$UpstreamControl_Phylop, alternative = "greater")

    Wilcoxon rank sum test with continuity correction

data:  NucReswPhyandC$phyloP and NucReswPhyandC$UpstreamControl_Phylop
W = 990510000, p-value < 2.2e-16
alternative hypothesis: true location shift is greater than 0

Actual are greater than region upstream

ggplot(NucReswPhyandCsmall, aes(x=SigPAU2, by=set, fill=set, y=Phylop)) + geom_boxplot() + scale_fill_brewer(palette = "Dark2",labels=c('PAS', 'Control') ) + stat_compare_means()

Version Author Date
7b73ce7 brimittleman 2020-03-04
NucReswPhyandCsmall_noc= NucReswPhyandCsmall %>% filter(set!="UpstreamControl_Phylop")

ggplot(NucReswPhyandCsmall_noc, aes(x=SigPAU2, fill=SigPAU2, y=Phylop )) + geom_boxplot() + stat_compare_means() + scale_fill_brewer(palette = "OrRd") + labs(title="Significant PAS are less conserved",x="")+ scale_x_discrete(labels=c("Not Significant", "Signficant"))+ theme(legend.position = "none",text= element_text(size=16))

Version Author Date
7b73ce7 brimittleman 2020-03-04

Significant are lower than not significant:

NucReswPhyandCsmall_nocYES= NucReswPhyandCsmall_noc %>% filter(SigPAU2=="Yes")
NucReswPhyandCsmall_nocNO= NucReswPhyandCsmall_noc %>% filter(SigPAU2=="No")

wilcox.test(NucReswPhyandCsmall_nocYES$Phylop, NucReswPhyandCsmall_nocNO$Phylop, alternative ="less")

    Wilcoxon rank sum test with continuity correction

data:  NucReswPhyandCsmall_nocYES$Phylop and NucReswPhyandCsmall_nocNO$Phylop
W = 55118000, p-value = 6.843e-06
alternative hypothesis: true location shift is less than 0

Significant have lower scores.

Negative enrichment

Number of negative in each set?

neg=NucReswPhyandCsmall %>% filter(Phylop <0) %>% group_by(set, SigPAU2) %>% summarise(nNeg=n())

pos=NucReswPhyandCsmall %>% filter(Phylop >0) %>% group_by(set, SigPAU2) %>% summarise(nPos=n())

both=neg %>% inner_join(pos,by= c('set', 'SigPAU2')) %>% mutate(PropNeg=nNeg/(nNeg+nPos))

both
# A tibble: 4 x 5
# Groups:   set [2]
  set                    SigPAU2  nNeg  nPos PropNeg
  <chr>                  <chr>   <int> <int>   <dbl>
1 phyloP                 No       8592 29079   0.228
2 phyloP                 Yes       788  2283   0.257
3 UpstreamControl_Phylop No      12085 25586   0.321
4 UpstreamControl_Phylop Yes      1062  2009   0.346

More negative overall in actual. Is there an enrichment for negative in the control set?

x=nrow(NucReswPhyandC %>% filter(SigPAU2=="Yes", UpstreamControl_Phylop<0))
m= nrow(NucReswPhyandC %>% filter(UpstreamControl_Phylop<0))
n=nrow(NucReswPhyandC %>% filter(UpstreamControl_Phylop>=0))
k=nrow(NucReswPhyandC %>% filter(SigPAU2=="Yes"))


#expected
which(grepl(max(dhyper(1:x, m, n, k)), dhyper(1:x, m, n, k)))
[1] 991
#actual:
x
[1] 1062
#pval
phyper(x,m,n,k,lower.tail=F)
[1] 0.002137693
b=nrow(NucReswPhyandC %>% filter(SigPAU2=="Yes", UpstreamControl_Phylop<0))
n=nrow(NucReswPhyandC %>% filter(SigPAU2=="Yes"))
B=nrow(NucReswPhyandC %>% filter(UpstreamControl_Phylop<0))
N=nrow(NucReswPhyandC)

(b/n)/(B/N)
[1] 1.071668

Stronger enrichement in the for negative in the real results compared to contol. 1.07x in control 1.11x in actual.

Maybe I need to move the control further up.

Is this a better control? Dont want to go into an exon? What about downstream?

Control sequence (downstream 200):

NucReswPhy_downstream=NucReswPhy %>% inner_join(metaStrand,by=c("chr", "start", "end")) %>% mutate(newStart=ifelse(strandFix=="+",  end, start-200), newEnd=ifelse(strandFix=="+", end+200, start))

NucReswPhy_downstreambed=NucReswPhy_downstream %>% select(chr, newStart, newEnd, PAS, Human, strandFix)

write.table(NucReswPhy_downstreambed,"../data/PhyloP/PAS_200downpregions.bed",col.names = F,row.names = F,quote = F,sep="\t")
python extractPhylopReg200down.py
Phylo200downCont=read.table("../data/PhyloP/PAS_phyloP_200downstream.txt",stringsAsFactors = F, col.names = c("chr", "start","end", "PAS","DownstreamControl_Phylop")) 

NucReswPhyandbothC=read.table("../data/DiffIso_Nuclear_DF/AllPAS_withGeneSig.txt", header = T, stringsAsFactors = F) %>% inner_join(phylores, by=c("chr", "start","end")) %>% inner_join(metaStrand,by=c("chr", "start", "end"))%>% inner_join(Phylo200UpContron, by="PAS")  %>% drop_na() %>%  inner_join(Phylo200downCont, by="PAS")  %>% drop_na()

NucReswPhyandCbothsmall=NucReswPhyandbothC %>% select(PAS,SigPAU2,phyloP ,UpstreamControl_Phylop,DownstreamControl_Phylop ) %>% gather("set", "Phylop", -PAS, -SigPAU2) %>% drop_na()


#difference in controls? 
wilcox.test(NucReswPhyandbothC$DownstreamControl_Phylop, NucReswPhyandbothC$UpstreamControl_Phylop,alternative = "greater")

    Wilcoxon rank sum test with continuity correction

data:  NucReswPhyandbothC$DownstreamControl_Phylop and NucReswPhyandbothC$UpstreamControl_Phylop
W = 771760000, p-value = 1
alternative hypothesis: true location shift is greater than 0
levels=NucReswPhyandCbothsmall$set %>% unique()
NucReswPhyandCbothsmall$set= factor(NucReswPhyandCbothsmall$set, levels = c("UpstreamControl_Phylop", "phyloP", "DownstreamControl_Phylop"))

my_comparisons <- list( c("DownsreamControl_Pylop", "phylopP"), c("DownsreamControl_Pylop", "UpstreamControl_Phylop"), c("phylopP", "UpstreamControl_Phylop") )
ggplot(NucReswPhyandCbothsmall, aes(x=set, by=set, fill=set, y=Phylop)) + geom_boxplot() + scale_fill_brewer(palette = "Dark2",labels=c("Upstream Control", "PAS", "Downstream Control") ) + stat_compare_means(ref.group = "phyloP",paired = FALSE,label = "p.signif")  + labs(x="", title="PAS conserved compared to surrounding regions" ) + scale_x_discrete( labels=c("Upstream Control", "PAS", "Downstream Control"))+ theme(legend.position = "none",text= element_text(size=16))

Version Author Date
0bbaafe brimittleman 2020-03-05

Same here. The actual region looks more conserved.

ggplot(NucReswPhyandCbothsmall, aes(x=SigPAU2, by=set, fill=set, y=Phylop)) + geom_boxplot() + scale_fill_brewer(palette = "Dark2",labels=c("Downstream Control", "PAS", "Upstream Control") )

Extend

2 more blocks up and downstream to add to plot.

Extend downstream:

NucReswPhy_downstream2=NucReswPhy_downstream %>% mutate(newStart2=ifelse(strandFix=="+",  newEnd, newStart-200), newEnd2=ifelse(strandFix=="+", newEnd+200, newStart))

NucReswPhy_downstream2bed=NucReswPhy_downstream2 %>% select(chr, newStart2, newEnd2, PAS, Human, strandFix)

write.table(NucReswPhy_downstream2bed,"../data/PhyloP/PAS_200downpregions2.bed",col.names = F,row.names = F,quote = F,sep="\t")

NucReswPhy_downstream3=NucReswPhy_downstream2 %>% mutate(newStart3=ifelse(strandFix=="+",  newEnd2, newStart2-200), newEnd3=ifelse(strandFix=="+", newEnd2+200, newStart2))

NucReswPhy_downstream3bed=NucReswPhy_downstream3 %>% select(chr, newStart3, newEnd3, PAS, Human, strandFix)

write.table(NucReswPhy_downstream3bed,"../data/PhyloP/PAS_200downpregions3.bed",col.names = F,row.names = F,quote = F,sep="\t")

Extend upstream:

NucReswPhy_upstream2=NucReswPhy_upstream %>%  mutate(newStart2=ifelse(strandFix=="+", newStart - 200, newEnd), newEnd2=ifelse(strandFix=="+", newStart, newEnd +200))

NucReswPhy_upstreambed2=NucReswPhy_upstream2 %>% select(chr, newStart2, newEnd2, PAS, Human, strandFix)

write.table(NucReswPhy_upstreambed2,"../data/PhyloP/PAS_200upregions2.bed",col.names = F,row.names = F,quote = F,sep="\t")

NucReswPhy_upstream3=NucReswPhy_upstream2 %>%  mutate(newStart3=ifelse(strandFix=="+", newStart2 - 200, newEnd2), newEnd3=ifelse(strandFix=="+", newStart2, newEnd2 +200))

NucReswPhy_upstreambed3=NucReswPhy_upstream3 %>% select(chr, newStart3, newEnd3, PAS, Human, strandFix)

write.table(NucReswPhy_upstreambed3,"../data/PhyloP/PAS_200upregions3.bed",col.names = F,row.names = F,quote = F,sep="\t")

Run phylop for each of these:

python extractPhylopGeneral.py ../data/PhyloP/PAS_200downpregions2.bed ../data/PhyloP/PAS_phyloP_200downstream2.txt

python extractPhylopGeneral.py ../data/PhyloP/PAS_200downpregions3.bed ../data/PhyloP/PAS_phyloP_200downstream3.txt

python extractPhylopGeneral.py ../data/PhyloP/PAS_200upregions2.bed ../data/PhyloP/PAS_phyloP_200upstream2.txt

python extractPhylopGeneral.py ../data/PhyloP/PAS_200upregions3.bed ../data/PhyloP/PAS_phyloP_200upstream3.txt
ResUpdown=NucReswPhyandbothC %>% select(PAS,SigPAU2,phyloP ,UpstreamControl_Phylop,DownstreamControl_Phylop ) 


Down2=read.table("../data/PhyloP/PAS_phyloP_200downstream2.txt",col.names = c("chr", "start", "end", "PAS", "Down2"),stringsAsFactors = F) %>% select(PAS, Down2)%>% drop_na()
Down3=read.table("../data/PhyloP/PAS_phyloP_200downstream3.txt",col.names = c("chr", "start", "end", "PAS", "Down3"),stringsAsFactors = F) %>% select(PAS, Down3)%>% drop_na()

Up2=read.table("../data/PhyloP/PAS_phyloP_200upstream2.txt",col.names = c("chr", "start", "end", "PAS", "Up2"),stringsAsFactors = F) %>% select(PAS, Up2)%>% drop_na()
Up3=read.table("../data/PhyloP/PAS_phyloP_200upstream3.txt",col.names = c("chr", "start", "end", "PAS", "Up3"),stringsAsFactors = F) %>% select(PAS, Up3)%>% drop_na()

ResUpdownAll= ResUpdown %>% inner_join(Down2, by="PAS")%>% inner_join(Down3, by="PAS") %>% inner_join(Up2, by="PAS") %>% inner_join(Up3, by="PAS")

ResUpdownAll_gather= ResUpdownAll %>% gather("Set", "PhyloP", -PAS, -SigPAU2)

ResUpdownAll_gather$Set=factor(ResUpdownAll_gather$Set, levels=c("Up3", "Up2","UpstreamControl_Phylop", "phyloP","DownstreamControl_Phylop", "Down2", "Down3" ))
ggplot(ResUpdownAll_gather, aes(x=Set, by=Set, fill=Set, y=PhyloP)) + geom_boxplot() + scale_fill_brewer(palette = "Dark2") + scale_x_discrete(labels=c("-600", "-400", "-200", '0','200','400','600')) + labs(x="Basepairs", title="PAS are more conserved than surrounding regions") + theme(legend.position = "none")

ggplot(ResUpdownAll_gather, aes(x=Set, by=Set, fill=Set, y=PhyloP)) + geom_boxplot() + scale_fill_brewer(palette = "Dark2") + facet_grid(~SigPAU2)+ scale_x_discrete(labels=c("-600", "-400", "-200", '0','200','400','600')) + labs(x="Basepairs", title="PAS are more conserved than surrounding regions") + theme(legend.position = "none")


sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] reshape2_1.4.3  ggpubr_0.2      magrittr_1.5    forcats_0.3.0  
 [5] stringr_1.3.1   dplyr_0.8.0.1   purrr_0.3.2     readr_1.3.1    
 [9] tidyr_0.8.3     tibble_2.1.1    ggplot2_3.1.1   tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5   haven_1.1.2        lattice_0.20-38   
 [4] colorspace_1.3-2   generics_0.0.2     htmltools_0.3.6   
 [7] yaml_2.2.0         utf8_1.1.4         rlang_0.4.0       
[10] later_0.7.5        pillar_1.3.1       glue_1.3.0        
[13] withr_2.1.2        RColorBrewer_1.1-2 modelr_0.1.2      
[16] readxl_1.1.0       plyr_1.8.4         munsell_0.5.0     
[19] gtable_0.2.0       workflowr_1.6.0    cellranger_1.1.0  
[22] rvest_0.3.2        evaluate_0.12      labeling_0.3      
[25] knitr_1.20         httpuv_1.4.5       fansi_0.4.0       
[28] broom_0.5.1        Rcpp_1.0.2         promises_1.0.1    
[31] scales_1.0.0       backports_1.1.2    jsonlite_1.6      
[34] fs_1.3.1           hms_0.4.2          digest_0.6.18     
[37] stringi_1.2.4      grid_3.5.1         rprojroot_1.3-2   
[40] cli_1.1.0          tools_3.5.1        lazyeval_0.2.1    
[43] crayon_1.3.4       whisker_0.3-2      pkgconfig_2.0.2   
[46] xml2_1.2.0         lubridate_1.7.4    assertthat_0.2.0  
[49] rmarkdown_1.10     httr_1.3.1         rstudioapi_0.10   
[52] R6_2.3.0           nlme_3.1-137       git2r_0.26.1      
[55] compiler_3.5.1