DEG_analysis_NTS_M_vs_F

workflowr

• Summary
• Checks
• Past versions

Last updated: 2023-12-23

Checks: 6 1

Knit directory: DO_Opioid/

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

---

R Markdown file: up-to-date

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.

Environment: empty

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.

Seed: set.seed(20200504)

The command set.seed(20200504) 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.

Session information: recorded

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

Cache: none

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

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute	relative
/projects/compsci/vmp/USERS/heh/DO_Opioid/data/BOT_NTS_rnaseq_results/	data/BOT_NTS_rnaseq_results

Repository version: 883216f

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 883216f. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

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:    .RData
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/Picture1.png

Untracked files:
    Untracked:  Rplot_rz.png
    Untracked:  TIMBR.test.prop.bm.MN.RET.RData
    Untracked:  TIMBR.test.random.RData
    Untracked:  TIMBR.test.rz.transformed_TVb_ml.RData
    Untracked:  analysis/DDO_morphine1_second_set_69k.stdout
    Untracked:  analysis/DO_Fentanyl.R
    Untracked:  analysis/DO_Fentanyl.err
    Untracked:  analysis/DO_Fentanyl.out
    Untracked:  analysis/DO_Fentanyl.sh
    Untracked:  analysis/DO_Fentanyl_69k.R
    Untracked:  analysis/DO_Fentanyl_69k.err
    Untracked:  analysis/DO_Fentanyl_69k.out
    Untracked:  analysis/DO_Fentanyl_69k.sh
    Untracked:  analysis/DO_Fentanyl_Cohort2_GCTA_herit.R
    Untracked:  analysis/DO_Fentanyl_Cohort2_gemma.R
    Untracked:  analysis/DO_Fentanyl_Cohort2_mapping.R
    Untracked:  analysis/DO_Fentanyl_Cohort2_mapping.err
    Untracked:  analysis/DO_Fentanyl_Cohort2_mapping.out
    Untracked:  analysis/DO_Fentanyl_Cohort2_mapping.sh
    Untracked:  analysis/DO_Fentanyl_GCTA_herit.R
    Untracked:  analysis/DO_Fentanyl_alternate_metrics_69k.R
    Untracked:  analysis/DO_Fentanyl_alternate_metrics_69k.err
    Untracked:  analysis/DO_Fentanyl_alternate_metrics_69k.out
    Untracked:  analysis/DO_Fentanyl_alternate_metrics_69k.sh
    Untracked:  analysis/DO_Fentanyl_alternate_metrics_array.R
    Untracked:  analysis/DO_Fentanyl_alternate_metrics_array.err
    Untracked:  analysis/DO_Fentanyl_alternate_metrics_array.out
    Untracked:  analysis/DO_Fentanyl_alternate_metrics_array.sh
    Untracked:  analysis/DO_Fentanyl_array.R
    Untracked:  analysis/DO_Fentanyl_array.err
    Untracked:  analysis/DO_Fentanyl_array.out
    Untracked:  analysis/DO_Fentanyl_array.sh
    Untracked:  analysis/DO_Fentanyl_combining2Cohort_GCTA_herit.R
    Untracked:  analysis/DO_Fentanyl_combining2Cohort_gemma.R
    Untracked:  analysis/DO_Fentanyl_combining2Cohort_mapping.R
    Untracked:  analysis/DO_Fentanyl_combining2Cohort_mapping.err
    Untracked:  analysis/DO_Fentanyl_combining2Cohort_mapping.out
    Untracked:  analysis/DO_Fentanyl_combining2Cohort_mapping.sh
    Untracked:  analysis/DO_Fentanyl_combining2Cohort_mapping_CoxPH.R
    Untracked:  analysis/DO_Fentanyl_finalreport_to_plink.sh
    Untracked:  analysis/DO_Fentanyl_gemma.R
    Untracked:  analysis/DO_Fentanyl_gemma.err
    Untracked:  analysis/DO_Fentanyl_gemma.out
    Untracked:  analysis/DO_Fentanyl_gemma.sh
    Untracked:  analysis/DO_morphine1.R
    Untracked:  analysis/DO_morphine1.Rout
    Untracked:  analysis/DO_morphine1.sh
    Untracked:  analysis/DO_morphine1.stderr
    Untracked:  analysis/DO_morphine1.stdout
    Untracked:  analysis/DO_morphine1_SNP.R
    Untracked:  analysis/DO_morphine1_SNP.Rout
    Untracked:  analysis/DO_morphine1_SNP.sh
    Untracked:  analysis/DO_morphine1_SNP.stderr
    Untracked:  analysis/DO_morphine1_SNP.stdout
    Untracked:  analysis/DO_morphine1_combined.R
    Untracked:  analysis/DO_morphine1_combined.Rout
    Untracked:  analysis/DO_morphine1_combined.sh
    Untracked:  analysis/DO_morphine1_combined.stderr
    Untracked:  analysis/DO_morphine1_combined.stdout
    Untracked:  analysis/DO_morphine1_combined_69k.R
    Untracked:  analysis/DO_morphine1_combined_69k.Rout
    Untracked:  analysis/DO_morphine1_combined_69k.sh
    Untracked:  analysis/DO_morphine1_combined_69k.stderr
    Untracked:  analysis/DO_morphine1_combined_69k.stdout
    Untracked:  analysis/DO_morphine1_combined_69k_m2.R
    Untracked:  analysis/DO_morphine1_combined_69k_m2.Rout
    Untracked:  analysis/DO_morphine1_combined_69k_m2.sh
    Untracked:  analysis/DO_morphine1_combined_69k_m2.stderr
    Untracked:  analysis/DO_morphine1_combined_69k_m2.stdout
    Untracked:  analysis/DO_morphine1_combined_weight_DOB.R
    Untracked:  analysis/DO_morphine1_combined_weight_DOB.Rout
    Untracked:  analysis/DO_morphine1_combined_weight_DOB.err
    Untracked:  analysis/DO_morphine1_combined_weight_DOB.out
    Untracked:  analysis/DO_morphine1_combined_weight_DOB.sh
    Untracked:  analysis/DO_morphine1_combined_weight_DOB.stderr
    Untracked:  analysis/DO_morphine1_combined_weight_DOB.stdout
    Untracked:  analysis/DO_morphine1_combined_weight_age.R
    Untracked:  analysis/DO_morphine1_combined_weight_age.err
    Untracked:  analysis/DO_morphine1_combined_weight_age.out
    Untracked:  analysis/DO_morphine1_combined_weight_age.sh
    Untracked:  analysis/DO_morphine1_combined_weight_age_GAMMT.R
    Untracked:  analysis/DO_morphine1_combined_weight_age_GAMMT.err
    Untracked:  analysis/DO_morphine1_combined_weight_age_GAMMT.out
    Untracked:  analysis/DO_morphine1_combined_weight_age_GAMMT.sh
    Untracked:  analysis/DO_morphine1_combined_weight_age_GAMMT_chr19.R
    Untracked:  analysis/DO_morphine1_combined_weight_age_GAMMT_chr19.err
    Untracked:  analysis/DO_morphine1_combined_weight_age_GAMMT_chr19.out
    Untracked:  analysis/DO_morphine1_combined_weight_age_GAMMT_chr19.sh
    Untracked:  analysis/DO_morphine1_cph.R
    Untracked:  analysis/DO_morphine1_cph.Rout
    Untracked:  analysis/DO_morphine1_cph.sh
    Untracked:  analysis/DO_morphine1_second_set.R
    Untracked:  analysis/DO_morphine1_second_set.Rout
    Untracked:  analysis/DO_morphine1_second_set.sh
    Untracked:  analysis/DO_morphine1_second_set.stderr
    Untracked:  analysis/DO_morphine1_second_set.stdout
    Untracked:  analysis/DO_morphine1_second_set_69k.R
    Untracked:  analysis/DO_morphine1_second_set_69k.Rout
    Untracked:  analysis/DO_morphine1_second_set_69k.sh
    Untracked:  analysis/DO_morphine1_second_set_69k.stderr
    Untracked:  analysis/DO_morphine1_second_set_SNP.R
    Untracked:  analysis/DO_morphine1_second_set_SNP.Rout
    Untracked:  analysis/DO_morphine1_second_set_SNP.sh
    Untracked:  analysis/DO_morphine1_second_set_SNP.stderr
    Untracked:  analysis/DO_morphine1_second_set_SNP.stdout
    Untracked:  analysis/DO_morphine1_second_set_weight_DOB.R
    Untracked:  analysis/DO_morphine1_second_set_weight_DOB.Rout
    Untracked:  analysis/DO_morphine1_second_set_weight_DOB.err
    Untracked:  analysis/DO_morphine1_second_set_weight_DOB.out
    Untracked:  analysis/DO_morphine1_second_set_weight_DOB.sh
    Untracked:  analysis/DO_morphine1_second_set_weight_DOB.stderr
    Untracked:  analysis/DO_morphine1_second_set_weight_DOB.stdout
    Untracked:  analysis/DO_morphine1_second_set_weight_age.R
    Untracked:  analysis/DO_morphine1_second_set_weight_age.Rout
    Untracked:  analysis/DO_morphine1_second_set_weight_age.err
    Untracked:  analysis/DO_morphine1_second_set_weight_age.out
    Untracked:  analysis/DO_morphine1_second_set_weight_age.sh
    Untracked:  analysis/DO_morphine1_second_set_weight_age.stderr
    Untracked:  analysis/DO_morphine1_second_set_weight_age.stdout
    Untracked:  analysis/DO_morphine1_weight_DOB.R
    Untracked:  analysis/DO_morphine1_weight_DOB.sh
    Untracked:  analysis/DO_morphine1_weight_age.R
    Untracked:  analysis/DO_morphine1_weight_age.sh
    Untracked:  analysis/DO_morphine_gemma.R
    Untracked:  analysis/DO_morphine_gemma.err
    Untracked:  analysis/DO_morphine_gemma.out
    Untracked:  analysis/DO_morphine_gemma.sh
    Untracked:  analysis/DO_morphine_gemma_firstmin.R
    Untracked:  analysis/DO_morphine_gemma_firstmin.err
    Untracked:  analysis/DO_morphine_gemma_firstmin.out
    Untracked:  analysis/DO_morphine_gemma_firstmin.sh
    Untracked:  analysis/DO_morphine_gemma_withpermu.R
    Untracked:  analysis/DO_morphine_gemma_withpermu.err
    Untracked:  analysis/DO_morphine_gemma_withpermu.out
    Untracked:  analysis/DO_morphine_gemma_withpermu.sh
    Untracked:  analysis/DO_morphine_gemma_withpermu_firstbatch_min.depression.R
    Untracked:  analysis/DO_morphine_gemma_withpermu_firstbatch_min.depression.err
    Untracked:  analysis/DO_morphine_gemma_withpermu_firstbatch_min.depression.out
    Untracked:  analysis/DO_morphine_gemma_withpermu_firstbatch_min.depression.sh
    Untracked:  analysis/Lisa_Tarantino_Interval_needs_mvar_annotation.R
    Untracked:  analysis/Plot_DO_morphine1_SNP.R
    Untracked:  analysis/Plot_DO_morphine1_SNP.Rout
    Untracked:  analysis/Plot_DO_morphine1_SNP.sh
    Untracked:  analysis/Plot_DO_morphine1_SNP.stderr
    Untracked:  analysis/Plot_DO_morphine1_SNP.stdout
    Untracked:  analysis/Plot_DO_morphine1_second_set_SNP.R
    Untracked:  analysis/Plot_DO_morphine1_second_set_SNP.Rout
    Untracked:  analysis/Plot_DO_morphine1_second_set_SNP.sh
    Untracked:  analysis/Plot_DO_morphine1_second_set_SNP.stderr
    Untracked:  analysis/Plot_DO_morphine1_second_set_SNP.stdout
    Untracked:  analysis/download_GSE100356_sra.sh
    Untracked:  analysis/fentanyl_2cohorts_coxph.R
    Untracked:  analysis/fentanyl_2cohorts_coxph.err
    Untracked:  analysis/fentanyl_2cohorts_coxph.out
    Untracked:  analysis/fentanyl_2cohorts_coxph.sh
    Untracked:  analysis/fentanyl_scanone.cph.R
    Untracked:  analysis/fentanyl_scanone.cph.err
    Untracked:  analysis/fentanyl_scanone.cph.out
    Untracked:  analysis/fentanyl_scanone.cph.sh
    Untracked:  analysis/geo_rnaseq.R
    Untracked:  analysis/heritability_first_second_batch.R
    Untracked:  analysis/morphine_fentanyl_survival_time.R
    Untracked:  analysis/nf-rnaseq-b6.R
    Untracked:  analysis/plot_fentanyl_2cohorts_coxph.R
    Untracked:  analysis/scripts/
    Untracked:  analysis/tibmr.R
    Untracked:  analysis/timbr_demo.R
    Untracked:  analysis/workflow_proc.R
    Untracked:  analysis/workflow_proc.sh
    Untracked:  analysis/workflow_proc.stderr
    Untracked:  analysis/workflow_proc.stdout
    Untracked:  analysis/x.R
    Untracked:  analysis/~.sh
    Untracked:  code/PLINKtoCSVR.R
    Untracked:  code/cfw/
    Untracked:  code/gemma_plot.R
    Untracked:  code/process.sanger.snp.R
    Untracked:  code/reconst_utils.R
    Untracked:  data/69k_grid_pgmap.RData
    Untracked:  data/BOT_NTS_rnaseq_results/
    Untracked:  data/CC_SARS-1/
    Untracked:  data/CC_SARS-2/
    Untracked:  data/Composite Post Kevins Program Group 2 Fentanyl Prepped for Hao.xlsx
    Untracked:  data/DO_WBP_Data_JAB_to_map.xlsx
    Untracked:  data/Fentanyl_alternate_metrics.xlsx
    Untracked:  data/FinalReport/
    Untracked:  data/GM/
    Untracked:  data/GM_covar.csv
    Untracked:  data/GM_covar_07092018_morphine.csv
    Untracked:  data/Jackson_Lab_Bubier_MURGIGV01/
    Untracked:  data/Lisa Tarantino Interval needs mvar.xlsx
    Untracked:  data/Lisa_Tarantino_Interval_needs_mvar_annotation.csv
    Untracked:  data/MPD_Upload_October.csv
    Untracked:  data/MPD_Upload_October_updated_sex.csv
    Untracked:  data/Master Fentanyl DO Study Sheet.xlsx
    Untracked:  data/MasterMorphine Second Set DO w DOB2.xlsx
    Untracked:  data/MasterMorphine Second Set DO.xlsx
    Untracked:  data/Morphine CC DO mice Updated with Published inbred strains.csv
    Untracked:  data/Morphine_CC_DO_mice_Updated_with_Published_inbred_strains.csv
    Untracked:  data/cc_variants.sqlite
    Untracked:  data/combined/
    Untracked:  data/fentanyl/
    Untracked:  data/fentanyl2/
    Untracked:  data/fentanyl_1_2/
    Untracked:  data/fentanyl_2cohorts_coxph_data.Rdata
    Untracked:  data/first/
    Untracked:  data/founder_geno.csv
    Untracked:  data/genetic_map.csv
    Untracked:  data/gm.json
    Untracked:  data/gwas.sh
    Untracked:  data/marker_grid_0.02cM_plus.txt
    Untracked:  data/metabolomics_mouse_fecal/
    Untracked:  data/mouse_genes_mgi.sqlite
    Untracked:  data/pheno.csv
    Untracked:  data/pheno_qtl2.csv
    Untracked:  data/pheno_qtl2_07092018_morphine.csv
    Untracked:  data/pheno_qtl2_w_dob.csv
    Untracked:  data/physical_map.csv
    Untracked:  data/rnaseq/
    Untracked:  data/sample_geno.csv
    Untracked:  data/second/
    Untracked:  figure/
    Untracked:  glimma-plots/
    Untracked:  head.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_MinDepressionRR_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_MinDepressionRR_coefplot_blup.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_RRDepressionRateHrSLOPE_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_RRRecoveryRateHrSLOPE_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_RRRecoveryRateHrSLOPE_coefplot_blup.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_StartofRecoveryHr_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_StartofRecoveryHr_coefplot_blup.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_Statusbin_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_Statusbin_coefplot_blup.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_SteadyStateDepressionDurationHrINTERVAL_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoDead(Hr)_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoDeadHr_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoDeadHr_coefplot_blup.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoProjectedRecoveryHr_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoProjectedRecoveryHr_coefplot_blup.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoSteadyRRDepression(Hr)_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoSteadyRRDepressionHr_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoSteadyRRDepressionHr_coefplot_blup.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoThresholdRecoveryHr_coefplot.pdf
    Untracked:  output/DO_Fentanyl_Cohort2_TimetoThresholdRecoveryHr_coefplot_blup.pdf
    Untracked:  output/DO_morphine_Min.depression.png
    Untracked:  output/DO_morphine_Min.depression22222_violin_chr5.pdf
    Untracked:  output/DO_morphine_Min.depression_coefplot.pdf
    Untracked:  output/DO_morphine_Min.depression_coefplot_blup.pdf
    Untracked:  output/DO_morphine_Min.depression_coefplot_blup_chr5.png
    Untracked:  output/DO_morphine_Min.depression_coefplot_blup_chrX.png
    Untracked:  output/DO_morphine_Min.depression_coefplot_chr5.png
    Untracked:  output/DO_morphine_Min.depression_coefplot_chrX.png
    Untracked:  output/DO_morphine_Min.depression_peak_genes_chr5.png
    Untracked:  output/DO_morphine_Min.depression_violin_chr5.png
    Untracked:  output/DO_morphine_Recovery.Time.png
    Untracked:  output/DO_morphine_Recovery.Time_coefplot.pdf
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_blup.pdf
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_blup_chr11.png
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_blup_chr4.png
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_blup_chr7.png
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_blup_chr9.png
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_chr11.png
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_chr4.png
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_chr7.png
    Untracked:  output/DO_morphine_Recovery.Time_coefplot_chr9.png
    Untracked:  output/DO_morphine_Status_bin.png
    Untracked:  output/DO_morphine_Status_bin_coefplot.pdf
    Untracked:  output/DO_morphine_Status_bin_coefplot_blup.pdf
    Untracked:  output/DO_morphine_Survival.Time.png
    Untracked:  output/DO_morphine_Survival.Time_coefplot.pdf
    Untracked:  output/DO_morphine_Survival.Time_coefplot_blup.pdf
    Untracked:  output/DO_morphine_Survival.Time_coefplot_blup_chr17.png
    Untracked:  output/DO_morphine_Survival.Time_coefplot_blup_chr8.png
    Untracked:  output/DO_morphine_Survival.Time_coefplot_chr17.png
    Untracked:  output/DO_morphine_Survival.Time_coefplot_chr8.png
    Untracked:  output/DO_morphine_combine_batch_peak_violin.pdf
    Untracked:  output/DO_morphine_combined_69k_m2_Min.depression.png
    Untracked:  output/DO_morphine_combined_69k_m2_Min.depression_coefplot.pdf
    Untracked:  output/DO_morphine_combined_69k_m2_Min.depression_coefplot_blup.pdf
    Untracked:  output/DO_morphine_combined_69k_m2_Recovery.Time.png
    Untracked:  output/DO_morphine_combined_69k_m2_Recovery.Time_coefplot.pdf
    Untracked:  output/DO_morphine_combined_69k_m2_Recovery.Time_coefplot_blup.pdf
    Untracked:  output/DO_morphine_combined_69k_m2_Status_bin.png
    Untracked:  output/DO_morphine_combined_69k_m2_Status_bin_coefplot.pdf
    Untracked:  output/DO_morphine_combined_69k_m2_Status_bin_coefplot_blup.pdf
    Untracked:  output/DO_morphine_combined_69k_m2_Survival.Time.png
    Untracked:  output/DO_morphine_combined_69k_m2_Survival.Time_coefplot.pdf
    Untracked:  output/DO_morphine_combined_69k_m2_Survival.Time_coefplot_blup.pdf
    Untracked:  output/DO_morphine_coxph_24hrs_kinship_QTL.png
    Untracked:  output/DO_morphine_cphout.RData
    Untracked:  output/DO_morphine_first_batch_peak_in_second_batch_violin.pdf
    Untracked:  output/DO_morphine_first_batch_peak_in_second_batch_violin_sidebyside.pdf
    Untracked:  output/DO_morphine_first_batch_peak_violin.pdf
    Untracked:  output/DO_morphine_operm.cph.RData
    Untracked:  output/DO_morphine_second_batch_on_first_batch_peak_violin.pdf
    Untracked:  output/DO_morphine_second_batch_peak_ch6surv_on_first_batchviolin.pdf
    Untracked:  output/DO_morphine_second_batch_peak_ch6surv_on_first_batchviolin2.pdf
    Untracked:  output/DO_morphine_second_batch_peak_in_first_batch_violin.pdf
    Untracked:  output/DO_morphine_second_batch_peak_in_first_batch_violin_sidebyside.pdf
    Untracked:  output/DO_morphine_second_batch_peak_violin.pdf
    Untracked:  output/DO_morphine_secondbatch_69k_Min.depression.png
    Untracked:  output/DO_morphine_secondbatch_69k_Min.depression_coefplot.pdf
    Untracked:  output/DO_morphine_secondbatch_69k_Min.depression_coefplot_blup.pdf
    Untracked:  output/DO_morphine_secondbatch_69k_Recovery.Time.png
    Untracked:  output/DO_morphine_secondbatch_69k_Recovery.Time_coefplot.pdf
    Untracked:  output/DO_morphine_secondbatch_69k_Recovery.Time_coefplot_blup.pdf
    Untracked:  output/DO_morphine_secondbatch_69k_Status_bin.png
    Untracked:  output/DO_morphine_secondbatch_69k_Status_bin_coefplot.pdf
    Untracked:  output/DO_morphine_secondbatch_69k_Status_bin_coefplot_blup.pdf
    Untracked:  output/DO_morphine_secondbatch_69k_Survival.Time.png
    Untracked:  output/DO_morphine_secondbatch_69k_Survival.Time_coefplot.pdf
    Untracked:  output/DO_morphine_secondbatch_69k_Survival.Time_coefplot_blup.pdf
    Untracked:  output/DO_morphine_secondbatch_Min.depression.png
    Untracked:  output/DO_morphine_secondbatch_Min.depression_coefplot.pdf
    Untracked:  output/DO_morphine_secondbatch_Min.depression_coefplot_blup.pdf
    Untracked:  output/DO_morphine_secondbatch_Recovery.Time.png
    Untracked:  output/DO_morphine_secondbatch_Recovery.Time_coefplot.pdf
    Untracked:  output/DO_morphine_secondbatch_Recovery.Time_coefplot_blup.pdf
    Untracked:  output/DO_morphine_secondbatch_Status_bin.png
    Untracked:  output/DO_morphine_secondbatch_Status_bin_coefplot.pdf
    Untracked:  output/DO_morphine_secondbatch_Status_bin_coefplot_blup.pdf
    Untracked:  output/DO_morphine_secondbatch_Survival.Time.png
    Untracked:  output/DO_morphine_secondbatch_Survival.Time_coefplot.pdf
    Untracked:  output/DO_morphine_secondbatch_Survival.Time_coefplot_blup.pdf
    Untracked:  output/Fentanyl/
    Untracked:  output/KPNA3.pdf
    Untracked:  output/SSC4D.pdf
    Untracked:  output/TIMBR.test.RData
    Untracked:  output/apr_69kchr_combined.RData
    Untracked:  output/apr_69kchr_k_loco_combined.rds
    Untracked:  output/apr_69kchr_second_set.RData
    Untracked:  output/combine_batch_variation.RData
    Untracked:  output/combined_gm.RData
    Untracked:  output/combined_gm.k_loco.rds
    Untracked:  output/combined_gm.k_overall.rds
    Untracked:  output/combined_gm.probs_8state.rds
    Untracked:  output/coxph/
    Untracked:  output/do.morphine.RData
    Untracked:  output/do.morphine.k_loco.rds
    Untracked:  output/do.morphine.probs_36state.rds
    Untracked:  output/do.morphine.probs_8state.rds
    Untracked:  output/do_Fentanyl_combine2cohort_MeanDepressionBR_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_MeanDepressionBR_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_MinDepressionBR_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_MinDepressionBR_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_MinDepressionRR_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_MinDepressionRR_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_RRRecoveryRateHr_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_RRRecoveryRateHr_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_StartofRecoveryHr_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_StartofRecoveryHr_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_Statusbin_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_Statusbin_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_SteadyStateDepressionDurationHr_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_SteadyStateDepressionDurationHr_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_SurvivalTime_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_SurvivalTime_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_TimetoDeadHr_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_TimetoDeadHr_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_TimetoMostlyDeadHr_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_TimetoMostlyDeadHr_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_TimetoProjectedRecoveryHr_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_TimetoProjectedRecoveryHr_coefplot_blup.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_TimetoRecoveryHr_coefplot.pdf
    Untracked:  output/do_Fentanyl_combine2cohort_TimetoRecoveryHr_coefplot_blup.pdf
    Untracked:  output/first_batch_variation.RData
    Untracked:  output/first_second_survival_peak_chr.xlsx
    Untracked:  output/hsq_1_first_batch_herit_qtl2.RData
    Untracked:  output/hsq_2_second_batch_herit_qtl2.RData
    Untracked:  output/morphine_fentanyl_survival_time.pdf
    Untracked:  output/old_temp/
    Untracked:  output/out_1_operm.RData
    Untracked:  output/pr_69kchr_combined.RData
    Untracked:  output/pr_69kchr_second_set.RData
    Untracked:  output/qtl.morphine.69k.out.combined.RData
    Untracked:  output/qtl.morphine.69k.out.combined_m2.RData
    Untracked:  output/qtl.morphine.69k.out.second_set.RData
    Untracked:  output/qtl.morphine.operm.RData
    Untracked:  output/qtl.morphine.out.RData
    Untracked:  output/qtl.morphine.out.combined_gm.RData
    Untracked:  output/qtl.morphine.out.combined_gm.female.RData
    Untracked:  output/qtl.morphine.out.combined_gm.male.RData
    Untracked:  output/qtl.morphine.out.combined_weight_DOB.RData
    Untracked:  output/qtl.morphine.out.combined_weight_age.RData
    Untracked:  output/qtl.morphine.out.female.RData
    Untracked:  output/qtl.morphine.out.male.RData
    Untracked:  output/qtl.morphine.out.second_set.RData
    Untracked:  output/qtl.morphine.out.second_set.female.RData
    Untracked:  output/qtl.morphine.out.second_set.male.RData
    Untracked:  output/qtl.morphine.out.second_set.weight_DOB.RData
    Untracked:  output/qtl.morphine.out.second_set.weight_age.RData
    Untracked:  output/qtl.morphine.out.weight_DOB.RData
    Untracked:  output/qtl.morphine.out.weight_age.RData
    Untracked:  output/qtl.morphine1.snpout.RData
    Untracked:  output/qtl.morphine2.snpout.RData
    Untracked:  output/sample_all_infor.csv
    Untracked:  output/sample_all_infor.txt
    Untracked:  output/second_batch_pheno.csv
    Untracked:  output/second_batch_variation.RData
    Untracked:  output/second_set_apr_69kchr_k_loco.rds
    Untracked:  output/second_set_gm.RData
    Untracked:  output/second_set_gm.k_loco.rds
    Untracked:  output/second_set_gm.probs_36state.rds
    Untracked:  output/second_set_gm.probs_8state.rds
    Untracked:  output/topSNP_chr5_mindepression.csv
    Untracked:  output/zoompeak_Min.depression_9.pdf
    Untracked:  output/zoompeak_Recovery.Time_16.pdf
    Untracked:  output/zoompeak_Status_bin_11.pdf
    Untracked:  output/zoompeak_Survival.Time_1.pdf
    Untracked:  output/zoompeak_fentanyl_Survival.Time_2.pdf
    Untracked:  sra-tools_v2.10.7.sif

Unstaged changes:
    Modified:   .gitignore
    Modified:   _workflowr.yml
    Deleted:    analysis/CC_SARS-2.Rmd
    Modified:   analysis/CC_SARS.Rmd
    Modified:   analysis/marker_violin.Rmd
    Modified:   output/CC_SARS_Chr16_QTL_interval.pdf
    Modified:   output/CC_SARS_Chr16_plotGeno.pdf
    Modified:   output/CC_SARS_Chr16_plotGeno.png

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 repository in which changes were made to the R Markdown (analysis/DEG_analysis_NTS_M_vs_F.Rmd) and HTML (docs/DEG_analysis_NTS_M_vs_F.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	883216f	xhyuo	2023-12-23	Sex DEG in tissue

---

workflow for rna-seq samples

library

library(stringr)
library(tidyverse)

Warning: replacing previous import 'lifecycle::last_warnings' by
'rlang::last_warnings' when loading 'hms'

Warning: replacing previous import 'ellipsis::check_dots_unnamed' by
'rlang::check_dots_unnamed' when loading 'hms'

Warning: replacing previous import 'ellipsis::check_dots_used' by
'rlang::check_dots_used' when loading 'hms'

Warning: replacing previous import 'ellipsis::check_dots_empty' by
'rlang::check_dots_empty' when loading 'hms'

Warning: package 'purrr' was built under R version 4.0.5

library(edgeR)
library(limma)
library(Glimma)
library(gplots)
library(org.Mm.eg.db)
library(RColorBrewer)
library(DESeq2)
library(pheatmap)
library(ggrepel)
library(DT)
library(enrichR)
library(cowplot)
library(ggplotify)
library(data.table)
library(WGCNA)
library(grid)
set.seed(123)

Collect RNA-seq emase out on gene level abundances.

#define the output directory
####merge all the genes.expected_read_counts####
all.dgerc.file <- list.files(path    = "/projects/compsci/vmp/USERS/heh/DO_Opioid/data/BOT_NTS_rnaseq_results/",
                             pattern = "\\.multiway.genes.expected_read_counts$", 
                             full.names = FALSE,
                             all.files  = TRUE,
                             recursive  = TRUE)
all.dgerc.file <- all.dgerc.file[!str_detect(all.dgerc.file, "Undetermined")]
#get the sample id
sampleid <- gsub(".*/", "", all.dgerc.file)
sampleid <- str_extract(sampleid, "^[^_]+_[^_]+")

#BOT 
BOT_sample <- map_dfr(
  c("/projects/activities/bubier/rnaseq/fastq/20230707_23-chesler-001/23-chesler-001_QCreport.csv","/projects/activities/bubier/rnaseq/fastq/20230710_23-chesler-001-run2/23-chesler-001-run2_QCreport.csv"), function(x){
    read_csv(x, skip = 16)
}) %>%
  dplyr::mutate(id = str_extract(Sample_Name, "^[^_]+_[^_]+"),
                .before = 1)

sheet <- readxl::read_excel("data/BOT_NTS_rnaseq_results/Bubier Brain Tracking Sheet.xlsx", sheet = 1, col_types = c("text", "text", "text", "numeric", "date", "logical", "text", "text", "text", "text", "text"), range = "A1:K393") %>%
  dplyr::select(1,2,3,5) %>%
  dplyr::distinct()

#data.frame
all.dgerc <- data.frame(file = all.dgerc.file, sampleid = sampleid) %>%
  dplyr::mutate(Tissue = case_when(
    (sampleid %in% BOT_sample$id) ~ "BOT",
    !(sampleid %in% BOT_sample$id) ~ "NTS"
  )) %>%
  dplyr::mutate(id = gsub("\\_GT23.*","",sampleid)) %>%
  dplyr::left_join(sheet, by = c("id" = "JCMS #"))

#39735
all.dgerc$Sex[all.dgerc$id == 39735] = c("F", "F")
all.dgerc$`Strain Name`[all.dgerc$id == 39735] = c("CC019/TauUncJ", "CC019/TauUncJ")
#40042
all.dgerc$Sex[all.dgerc$id == 40042] = c("M", "M")
all.dgerc$`Strain Name`[all.dgerc$id == 40042] = c("CC013/GeniUncJ", "CC013/GeniUncJ")
#40043
all.dgerc$Sex[all.dgerc$id == 40043] = "M"
all.dgerc$`Strain Name`[all.dgerc$id == 40043] = "CC013/GeniUncJ"

#merge
command.merge.dgerc <- paste0("bash -c 'paste ", paste(paste0("<(cut -f 10 /projects/compsci/vmp/USERS/heh/DO_Opioid/data/BOT_NTS_rnaseq_results/", all.dgerc.file,")"), collapse = " "), " > /projects/compsci/vmp/USERS/heh/DO_Opioid/data/BOT_NTS_rnaseq_results/merged.dgerc'")
system(command.merge.dgerc)

#read into R
merged.dgerc <- read.table("data/BOT_NTS_rnaseq_results/merged.dgerc", header = T, sep = "\t")
colnames(merged.dgerc) <- sampleid
expgene <- read.table(file = paste0("data/BOT_NTS_rnaseq_results/",as.character(all.dgerc$file[1])),header=TRUE,sep="\t")
rownames(merged.dgerc) <- expgene[,1]
write.csv(merged.dgerc,file="data/BOT_NTS_rnaseq_results/merged.dgerc.csv",quote=F,row.names=T)

save(merged.dgerc, file="data/BOT_NTS_rnaseq_results/merged.dgerc.RData")
save(all.dgerc, file="data/BOT_NTS_rnaseq_results/all.dgerc.RData")

count matrix and design matrix

load("data/BOT_NTS_rnaseq_results/all.dgerc.RData")
#RNA seq count data
countdata <- get(load("data/BOT_NTS_rnaseq_results/merged.dgerc.RData"))
#floor countdata
countdata <- floor(countdata)
#Removing genes that are lowly expressed as 0
countdata <- countdata[rowSums(countdata) != 0,]

#gene annotation
genes <- AnnotationDbi::select(org.Mm.eg.db, keys=rownames(countdata), columns=c("SYMBOL"), 
                               keytype="ENSEMBL")

'select()' returned 1:many mapping between keys and columns

genes <- genes[!duplicated(genes$ENSEMBL),]
genes <- genes[match(rownames(countdata), genes$ENSEMBL),]
#order
all.equal(rownames(countdata), genes$ENSEMBL)

[1] TRUE

#"35714_NTS_M" "38604_NTS_M" "38663_NTS_M" "40073_NTS_M"
#design matrix
design.matrix <- all.dgerc %>%
  dplyr::select(-1) %>%
  dplyr::rename(Strain = 5) %>%
  dplyr::mutate(name = case_when(
    id %in% c(35714, 38604, 38663, 40073) ~ paste(all.dgerc$sampleid, all.dgerc$Tissue, all.dgerc$Sex, sep = "_"),
    TRUE ~ paste(all.dgerc$id, all.dgerc$Tissue, all.dgerc$Sex, sep = "_")
  )) %>%
  group_by(name) %>% 
  dplyr::mutate(name = make.unique(as.character(name))) %>%
  dplyr::mutate(across(c(Sex, Strain, Tissue), as.factor))
rownames(design.matrix) = design.matrix$name

Warning: Setting row names on a tibble is deprecated.

#strain_means
strain_means = read_csv("data/BOT_NTS_rnaseq_results/Project1046/Project1046_strainmeans.csv")

── Column specification ────────────────────────────────────────────────────────
cols(
  measnum = col_double(),
  varname = col_character(),
  strain = col_character(),
  strainid = col_double(),
  sex = col_character(),
  mean = col_double(),
  nmice = col_double(),
  sd = col_double(),
  sem = col_double(),
  cv = col_double(),
  minval = col_double(),
  maxval = col_double(),
  zscore = col_double()
)

strain_means_wide = strain_means %>%
  dplyr::select(2,3,5,6) %>%
  dplyr::filter(str_detect(varname, "Baseline|Time|Depression|Recovery")) %>%
  tidyr::pivot_wider(names_from = varname, values_from = mean) %>%
  dplyr::rename(Strain = strain, Sex = sex) %>%
  dplyr::mutate(Sex = toupper(Sex))
#left_join
design.matrix.pheno = design.matrix %>%
  left_join(strain_means_wide) %>%
  dplyr::select(where(~ !all(is.na(.)))) %>%
  as.data.frame()

Joining, by = c("Sex", "Strain")

rownames(design.matrix.pheno) = design.matrix.pheno$name


#design.matrix$group <- factor(paste0(design.matrix$Strain, design.matrix$Tissue))
#order
all.equal(design.matrix$sampleid, colnames(countdata))

[1] TRUE

#new colname of countdata
colnames(countdata) = rownames(design.matrix)

#subset to BOT
design.matrix <- design.matrix[design.matrix$Tissue == "NTS",]
rownames(design.matrix) = design.matrix$name

Warning: Setting row names on a tibble is deprecated.

countdata <- countdata[, rownames(design.matrix)]
all.equal(colnames(countdata), rownames(design.matrix))

[1] TRUE

#To now construct the DESeqDataSet object from the matrix of counts and the sample information table, we use:
ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
                                 colData   = design.matrix,
                                 design    = ~Strain + Sex)

converting counts to integer mode

  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not a warning or an error]

QC process

#Pre-filtering the dataset
#perform a minimal pre-filtering to keep only rows that have at least 10 reads total.
keep <- rowSums(counts(ddsMat) >= 10) >= 2
ddsMat <- ddsMat[keep,]
ddsMat

class: DESeqDataSet 
dim: 21148 188 
metadata(1): version
assays(1): counts
rownames(21148): ENSMUSG00000000001 ENSMUSG00000000028 ...
  ENSMUSG00001074846 ENSMUSG00002076083
rowData names(0):
colnames(188): 34303_NTS_F 34305_NTS_M ... 41413_NTS_F 41415_NTS_F
colData names(7): sampleid Tissue ... DOB name

nrow(ddsMat)

[1] 21148

## [1] 27095

# DESeq2 creates a matrix when you use the counts() function
## First convert normalized_counts to a data frame and transfer the row names to a new column called "gene"
# this gives log2(n + 1)
ntd <- normTransform(ddsMat)

  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not a warning or an error]

normalized_counts <- assay(ntd) %>% 
  data.frame() %>%
  rownames_to_column(var="gene") %>%
  as_tibble()

#The variance stabilizing transformation and the rlog
#The rlog tends to work well on small datasets (n < 30), potentially outperforming the VST when there is a wide range of sequencing depth across samples (an order of magnitude difference). 
#rld <- vst(ddsMat, blind = FALSE)
#save(rld, file = "data/BOT_NTS_rnaseq_results/NTS_rld.RData")
load("data/BOT_NTS_rnaseq_results/NTS_rld.RData")
#head(assay(rld), 3)
#sample distance
sampleDists <- dist(t(assay(rld)))
sampleDistMatrix <- as.matrix( sampleDists )
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
#annotation
df <- as.data.frame(colData(ddsMat)[, c("Strain", "Sex")])
#heatmap on sample distance
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors,
         annotation_col  = df, 
         show_rownames = F, show_colnames = F,
         annotation_colors = list(
           Sex = c("M" = "#0064C9", 
                   "F"    ="#B10DC9")),
         border_color = NA)

#heatmap on correlation matrix
rld_cor <- cor(assay(rld))
pheatmap(rld_cor,
         annotation_col  = df, 
         clustering_distance_rows = "correlation",
         clustering_distance_cols = "correlation",
         show_rownames = F, show_colnames = F,
         annotation_colors = list(
                                  Sex = c("M" = "#0064C9", 
                                          "F"    ="#B10DC9")),
         border_color = NA)

#PCA plot
#Another way to visualize sample-to-sample distances is a principal components analysis (PCA). 
pca.plot <- plotPCA(rld, intgroup = c("Sex"), returnData = FALSE)
pca.plot

Differential expression analysis without interaction term

design(ddsMat)

~Strain + Sex

#Running the differential expression pipeline
#res <- DESeq(ddsMat)
#save(res, file = "data/BOT_NTS_rnaseq_results/NTS_res_no_interaction.RData")
load("data/BOT_NTS_rnaseq_results/NTS_res_no_interaction.RData")
resultsNames(res)

 [1] "Intercept"                              
 [2] "Strain_CC059.TauUncJ_vs_CC023.GeniUncJ" 
 [3] "Strain_CC030.GeniUncJ_vs_CC023.GeniUncJ"
 [4] "Strain_CC079.TauUncJ_vs_CC023.GeniUncJ" 
 [5] "Strain_CC010.GeniUncJ_vs_CC023.GeniUncJ"
 [6] "Strain_CC009.UncJ_vs_CC023.GeniUncJ"    
 [7] "Strain_CC011.UncJ_vs_CC023.GeniUncJ"    
 [8] "Strain_CC016.GeniUncJ_vs_CC023.GeniUncJ"
 [9] "Strain_CC024.GeniUncJ_vs_CC023.GeniUncJ"
[10] "Strain_CC065.UncJ_vs_CC023.GeniUncJ"    
[11] "Strain_CC043.GeniUncJ_vs_CC023.GeniUncJ"
[12] "Strain_CC019.TauUncJ_vs_CC023.GeniUncJ" 
[13] "Strain_CC080.TauUncJ_vs_CC023.GeniUncJ" 
[14] "Strain_CC084.TauJ_vs_CC023.GeniUncJ"    
[15] "Strain_CC027.GeniUncJ_vs_CC023.GeniUncJ"
[16] "Strain_CC028.GeniUncJ_vs_CC023.GeniUncJ"
[17] "Strain_CC004.TauUncJ_vs_CC023.GeniUncJ" 
[18] "Strain_CC002.UncJ_vs_CC023.GeniUncJ"    
[19] "Strain_CC008.GeniUncJ_vs_CC023.GeniUncJ"
[20] "Strain_CC051.TauUncJ_vs_CC023.GeniUncJ" 
[21] "Strain_CC020.GeniUncJ_vs_CC023.GeniUncJ"
[22] "Strain_CC032.GeniUncJ_vs_CC023.GeniUncJ"
[23] "Strain_CC001.UncJ_vs_CC023.GeniUncJ"    
[24] "Strain_CC078.TauUncJ_vs_CC023.GeniUncJ" 
[25] "Strain_CC061.GeniUncJ_vs_CC023.GeniUncJ"
[26] "Strain_CC041.TauUncJ_vs_CC023.GeniUncJ" 
[27] "Strain_CC029.UncJ_vs_CC023.GeniUncJ"    
[28] "Strain_CC018.UncJ_vs_CC023.GeniUncJ"    
[29] "Strain_CC040.TauUncJ_vs_CC023.GeniUncJ" 
[30] "Strain_CC037.TauUncJ_vs_CC023.GeniUncJ" 
[31] "Strain_CC060.UncJ_vs_CC023.GeniUncJ"    
[32] "Strain_CC033.GeniUncJ_vs_CC023.GeniUncJ"
[33] "Strain_CC074.UncJ_vs_CC023.GeniUncJ"    
[34] "Strain_CC035.UncJ_vs_CC023.GeniUncJ"    
[35] "Strain_CC058.UncJ_vs_CC023.GeniUncJ"    
[36] "Strain_CC044.UncJ_vs_CC023.GeniUncJ"    
[37] "Strain_CC015.UncJ_vs_CC023.GeniUncJ"    
[38] "Strain_CC036.UncJ_vs_CC023.GeniUncJ"    
[39] "Strain_CC012.GeniUncJ_vs_CC023.GeniUncJ"
[40] "Strain_CC013.GeniUncJ_vs_CC023.GeniUncJ"
[41] "Strain_CC026.GeniUncJ_vs_CC023.GeniUncJ"
[42] "Strain_CC083.UncJ_vs_CC023.GeniUncJ"    
[43] "Strain_CC003.UncJ_vs_CC023.GeniUncJ"    
[44] "Strain_CC017.UncJ_vs_CC023.GeniUncJ"    
[45] "Strain_CC025.GeniUncJ_vs_CC023.GeniUncJ"
[46] "Strain_CC057.UncJ_vs_CC023.GeniUncJ"    
[47] "Strain_CC005.TauUncJ_vs_CC023.GeniUncJ" 
[48] "Strain_CC007.UncJ_vs_CC023.GeniUncJ"    
[49] "Strain_CC068.TauUncJ_vs_CC023.GeniUncJ" 
[50] "Strain_CC045.GeniUncJ_vs_CC023.GeniUncJ"
[51] "Strain_CC075.UncJ_vs_CC023.GeniUncJ"    
[52] "Sex_M_vs_F"                             

#comparison for Sex------
#Building the results table
res.tab.sex <- results(res, name = "Sex_M_vs_F", alpha = 0.05)
res.tab.sex

log2 fold change (MLE): Sex M vs F 
Wald test p-value: Sex M vs F 
DataFrame with 21148 rows and 6 columns
                    baseMean log2FoldChange     lfcSE       stat    pvalue
                   <numeric>      <numeric> <numeric>  <numeric> <numeric>
ENSMUSG00000000001 976.40554     0.00234825 0.0249332  0.0941819  0.924965
ENSMUSG00000000028  36.23281     0.04105280 0.0540986  0.7588513  0.447942
ENSMUSG00000000031  12.14730     0.12133660 0.2141039  0.5667183  0.570906
ENSMUSG00000000037  39.92724    -0.05348666 0.0700515 -0.7635332  0.445145
ENSMUSG00000000049   6.00797     0.03056583 0.1415496  0.2159372  0.829037
...                      ...            ...       ...        ...       ...
ENSMUSG00000118669 789.57162     -0.0322927 0.0311901  -1.035353 0.3005041
ENSMUSG00000118670   7.42688      0.1738268 0.3075212   0.565251 0.5719027
ENSMUSG00000118671  43.01866      0.2766873 0.1582934   1.747940 0.0804744
ENSMUSG00001074846   1.12849     -0.0905424 0.4217898  -0.214662 0.8300306
ENSMUSG00002076083 861.86663      0.0111824 0.0223017   0.501414 0.6160796
                        padj
                   <numeric>
ENSMUSG00000000001  0.999799
ENSMUSG00000000028  0.977535
ENSMUSG00000000031  0.991496
ENSMUSG00000000037  0.976891
ENSMUSG00000000049  0.998322
...                      ...
ENSMUSG00000118669  0.965828
ENSMUSG00000118670  0.991496
ENSMUSG00000118671  0.873792
ENSMUSG00001074846  0.998322
ENSMUSG00002076083  0.991496

summary(res.tab.sex)

out of 21148 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 10, 0.047%
LFC < 0 (down)     : 17, 0.08%
outliers [1]       : 45, 0.21%
low counts [2]     : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

table(res.tab.sex$padj < 0.05)

FALSE  TRUE 
21076    27 

#We subset the results table to these genes and then sort it by the log2 fold change estimate to get the significant genes with the strongest down-regulation:
resSig.sex <- subset(res.tab.sex, padj < 0.05)
head(resSig.sex[order(resSig.sex$log2FoldChange), ])

log2 fold change (MLE): Sex M vs F 
Wald test p-value: Sex M vs F 
DataFrame with 6 rows and 6 columns
                     baseMean log2FoldChange     lfcSE      stat      pvalue
                    <numeric>      <numeric> <numeric> <numeric>   <numeric>
ENSMUSG00000085715   58.51883      -4.080040 0.2277426 -17.91513 8.98519e-72
ENSMUSG00000048173    5.80514      -2.445916 0.5865094  -4.17029 3.04209e-05
ENSMUSG00000040533    4.92650      -0.664170 0.1656687  -4.00903 6.09695e-05
ENSMUSG00000037369 1138.57411      -0.527971 0.0301824 -17.49269 1.62883e-68
ENSMUSG00000035150 6066.43743      -0.511658 0.0241404 -21.19509 1.05999e-99
ENSMUSG00000020241  187.15543      -0.459576 0.0787469  -5.83612 5.34316e-09
                          padj
                     <numeric>
ENSMUSG00000085715 2.70878e-68
ENSMUSG00000048173 2.56789e-02
ENSMUSG00000040533 4.76533e-02
ENSMUSG00000037369 4.29665e-65
ENSMUSG00000035150 4.47378e-96
ENSMUSG00000020241 7.51712e-06

# with the strongest up-regulation:
head(resSig.sex[order(resSig.sex$log2FoldChange, decreasing = TRUE), ])

log2 fold change (MLE): Sex M vs F 
Wald test p-value: Sex M vs F 
DataFrame with 6 rows and 6 columns
                     baseMean log2FoldChange     lfcSE      stat       pvalue
                    <numeric>      <numeric> <numeric> <numeric>    <numeric>
ENSMUSG00000069045 1103.48013       15.44290  0.475455  32.48024 2.02735e-231
ENSMUSG00000069044    4.65902       15.02171  2.995862   5.01415  5.32681e-07
ENSMUSG00000069049  406.65733       13.56559  0.448186  30.26777 3.04591e-201
ENSMUSG00000068457  512.20277       11.93478  0.411961  28.97061 1.54362e-184
ENSMUSG00000056673  613.08307        9.98412  0.355190  28.10924 7.55317e-174
ENSMUSG00000099876   17.46157        7.00894  0.357500  19.60542  1.39002e-85
                           padj
                      <numeric>
ENSMUSG00000069045 4.27832e-227
ENSMUSG00000069044  6.61245e-04
ENSMUSG00000069049 3.21390e-197
ENSMUSG00000068457 1.08584e-180
ENSMUSG00000056673 3.98486e-170
ENSMUSG00000099876  4.88892e-82

Visualization

#Visualization for Sex result------
#Volcano plot
## Obtain logical vector regarding whether padj values are less than 0.05
threshold_OE <- (res.tab.sex$padj < 0.05 & abs(res.tab.sex$log2FoldChange) >= 1)
## Determine the number of TRUE values
length(which(threshold_OE))

[1] 9

## Add logical vector as a column (threshold) to the res.tab.sex
res.tab.sex$threshold <- threshold_OE 
## Sort by ordered padj
res.tab.sex_ordered <- res.tab.sex %>%
  data.frame() %>%
  rownames_to_column(var="ENSEMBL") %>% 
  arrange(padj) %>%
  mutate(genelabels = "") %>%
  as_tibble() %>%
  left_join(genes)

Joining, by = "ENSEMBL"

## Create a column to indicate which genes to label
res.tab.sex_ordered$genelabels[1:10] <- res.tab.sex_ordered$SYMBOL[1:10]

#display res.tab.sex_ordered
DT::datatable(res.tab.sex_ordered[res.tab.sex_ordered$padj < 0.05,],
              filter = list(position = 'top', clear = FALSE),
              extensions = 'Buttons',
              options = list(dom = 'Blfrtip',
                         buttons = c('csv', 'excel'),
                         lengthMenu = list(c(10,25,50,-1),
                                           c(10,25,50,"All")),
                         pageLength = 40, 
                             scrollY = "300px", 
                             scrollX = "40px"),
              caption = htmltools::tags$caption(style = 'caption-side: top; text-align: left; color:black; font-size:200% ;','DEG analysis sex results (fdr < 0.05)'))

#Volcano plot
volcano.plot.sex <- ggplot(res.tab.sex_ordered) +
  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = threshold)) +
  scale_color_manual(values=c("blue", "red")) +
  geom_text_repel(aes(x = log2FoldChange, y = -log10(padj), 
                      label = genelabels,
                      size = 3.5)) +
  ggtitle("DEG for sex Male vs Female") +
  xlab("log2 fold change") + 
  ylab("-log10 adjusted p-value") +
  theme(legend.position = "none",
        plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = rel(1.25))) 
#print(volcano.plot.sex)

#heatmap
# Extract normalized expression for significant genes fdr < 0.05 & abs(log2FoldChange) >= 1)
normalized_counts_sig.sex <- normalized_counts %>% 
  filter(gene %in% rownames(subset(resSig.sex, padj < 0.05 & abs(log2FoldChange) >= 1)))

### Set a color palette
heat_colors <- brewer.pal(6, "YlOrRd")
#annotation
df <- as.data.frame(colData(ddsMat)[,c("Strain", "Sex")])
### Run pheatmap using the metadata data frame for the annotation
dat = as.matrix(normalized_counts_sig.sex[,-1])
colnames(dat) = rownames(df)
pheatmap(dat, 
         color = heat_colors, 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col  = df, 
         annotation_colors = list(
                                  Sex = c("M" = "#0064C9", 
                                          "F"    ="#B10DC9")),
         border_color = NA, 
         fontsize = 10, 
         scale = "row", 
         fontsize_row = 10, 
         height = 20,
         main = "Heatmap of the top DEGs in Sex Male vs Femal(fdr < 0.05 & abs(log2FoldChange) >= 1)")

Enrichment

#enrichment analysis
dbs <- c("WikiPathways_2019_Mouse",
         "GO_Biological_Process_2021",
         "GO_Cellular_Component_2021",
         "GO_Molecular_Function_2021",
         "KEGG_2019_Mouse",
         "Mouse_Gene_Atlas",
         "MGI_Mammalian_Phenotype_Level_4_2019")

#sex results------
resSig.sex.tab <- resSig.sex %>%
  data.frame() %>%
  rownames_to_column(var="ENSEMBL") %>%
  as_tibble() %>%
  left_join(genes)

Joining, by = "ENSEMBL"

#up-regulated sex genes---------------------------
up.genes <- resSig.sex.tab %>%
  filter(log2FoldChange > 0) %>%
  pull(SYMBOL)
#up-regulated genes enrichment
up.genes.enriched <- enrichr(as.character(na.omit(up.genes)), dbs)

Uploading data to Enrichr... Done.
  Querying WikiPathways_2019_Mouse... Done.
  Querying GO_Biological_Process_2021... Done.
  Querying GO_Cellular_Component_2021... Done.
  Querying GO_Molecular_Function_2021... Done.
  Querying KEGG_2019_Mouse... Done.
  Querying Mouse_Gene_Atlas... Done.
  Querying MGI_Mammalian_Phenotype_Level_4_2019... Done.
Parsing results... Done.

for (j in 1:length(up.genes.enriched)){
  up.genes.enriched[[j]] <- cbind(data.frame(Library = names(up.genes.enriched)[j]),up.genes.enriched[[j]])
}
up.genes.enriched <- do.call(rbind.data.frame, up.genes.enriched) %>%
  filter(Adjusted.P.value <= 0.05) %>%
  mutate(logpvalue = -log10(P.value)) %>%
  arrange(desc(logpvalue))
#display up.genes.enriched
DT::datatable(up.genes.enriched,filter = list(position = 'top', clear = FALSE),
              extensions = 'Buttons',
              options = list(dom = 'Blfrtip',
                         buttons = c('csv', 'excel'),
                         lengthMenu = list(c(10,25,50,-1),
                                           c(10,25,50,"All")),
                         pageLength = 40, 
                             scrollY = "300px", 
                             scrollX = "40px")
              )

#barpot
up.genes.enriched.plot <- up.genes.enriched %>%
  filter(Adjusted.P.value <= 0.05) %>%
  mutate(Term = fct_reorder(Term, -logpvalue)) %>%
  ggplot(data = ., aes(x = Term, y = logpvalue, fill = Library, label = Overlap)) +
  geom_bar(stat = "identity") +
  geom_text(position = position_dodge(width = 0.9),
            hjust = 0) + 
  theme_bw() + 
  ylab("-log(p.value)") +
  xlab("") + 
  ggtitle("Enrichment of up-regulated sex genes \n(Adjusted.P.value <= 0.05)") +
  theme(plot.background = element_blank() ,
        panel.border = element_blank(), 
        panel.background = element_blank(),
        #legend.position = "none",
        plot.title = element_text(hjust = 0),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) + 
  theme(axis.line = element_line(color = 'black')) + 
  theme(axis.title.x = element_text(size = 12, vjust=-0.5)) + 
  theme(axis.title.y = element_text(size = 12, vjust= 1.0)) + 
  theme(axis.text = element_text(size = 12)) + 
  theme(plot.title = element_text(size = 12)) + 
  coord_flip()
up.genes.enriched.plot

#down-regulated sex genes---------------------------
down.genes <- resSig.sex.tab %>%
  filter(log2FoldChange < 0) %>%
  pull(SYMBOL)
#down-regulated genes enrichment
down.genes.enriched <- enrichr(as.character(na.omit(down.genes)), dbs)

Uploading data to Enrichr... Done.
  Querying WikiPathways_2019_Mouse... Done.
  Querying GO_Biological_Process_2021... Done.
  Querying GO_Cellular_Component_2021... Done.
  Querying GO_Molecular_Function_2021... Done.
  Querying KEGG_2019_Mouse... Done.
  Querying Mouse_Gene_Atlas... Done.
  Querying MGI_Mammalian_Phenotype_Level_4_2019... Done.
Parsing results... Done.

for (j in 1:length(down.genes.enriched)){
  down.genes.enriched[[j]] <- cbind(data.frame(Library = names(down.genes.enriched)[j]),down.genes.enriched[[j]])
}
down.genes.enriched <- do.call(rbind.data.frame, down.genes.enriched) %>%
  filter(Adjusted.P.value <= 0.05) %>%
  mutate(logpvalue = -log10(P.value)) %>%
  arrange(desc(logpvalue)) 
#display down.genes.enriched
DT::datatable(down.genes.enriched,
              filter = list(position = 'top', clear = FALSE),
              extensions = 'Buttons',
              options = list(dom = 'Blfrtip',
                         buttons = c('csv', 'excel'),
                         lengthMenu = list(c(10,25,50,-1),
                                           c(10,25,50,"All")),
                         pageLength = 40, 
                             scrollY = "300px", 
                             scrollX = "40px")
              )

#barpot
down.genes.enriched.plot <- down.genes.enriched %>%
  filter(Adjusted.P.value <= 0.05) %>%
  mutate(Term = fct_reorder(Term, -logpvalue)) %>%
  ggplot(data = ., aes(x = Term, y = logpvalue, fill = Library, label = Overlap)) +
  geom_bar(stat = "identity") +
  geom_text(position = position_dodge(width = 0.9),
            hjust = 0) + 
  theme_bw() + 
  ylab("-log(p.value)") +
  xlab("") + 
  ggtitle("Enrichment of down-regulated sex genes \n(Adjusted.P.value <= 0.05)") +
  theme(plot.background = element_blank() ,
        panel.border = element_blank(), 
        panel.background = element_blank(),
        #legend.position = "none",
        plot.title = element_text(hjust = 0),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) + 
  theme(axis.line = element_line(color = 'black')) + 
  theme(axis.title.x = element_text(size = 12, vjust=-0.5)) + 
  theme(axis.title.y = element_text(size = 12, vjust= 1.0)) + 
  theme(axis.text = element_text(size = 12)) + 
  theme(plot.title = element_text(size = 12)) + 
  coord_flip()
down.genes.enriched.plot

Session information

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.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=C             
 [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] grid      parallel  stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] WGCNA_1.69                  fastcluster_1.1.25         
 [3] dynamicTreeCut_1.63-1       data.table_1.13.6          
 [5] ggplotify_0.0.5             cowplot_1.1.1              
 [7] enrichR_3.0                 DT_0.17                    
 [9] ggrepel_0.9.1               pheatmap_1.0.12            
[11] DESeq2_1.30.1               SummarizedExperiment_1.20.0
[13] MatrixGenerics_1.2.1        matrixStats_0.58.0         
[15] GenomicRanges_1.42.0        GenomeInfoDb_1.26.7        
[17] RColorBrewer_1.1-2          org.Mm.eg.db_3.12.0        
[19] AnnotationDbi_1.52.0        IRanges_2.24.1             
[21] S4Vectors_0.28.1            Biobase_2.50.0             
[23] BiocGenerics_0.36.1         gplots_3.1.1               
[25] Glimma_2.0.0                edgeR_3.32.1               
[27] limma_3.46.0                forcats_0.5.1              
[29] dplyr_1.0.4                 purrr_0.3.4                
[31] readr_1.4.0                 tidyr_1.1.2                
[33] tibble_3.0.6                ggplot2_3.3.3              
[35] tidyverse_1.3.0             stringr_1.4.0              
[37] workflowr_1.6.2            

loaded via a namespace (and not attached):
  [1] readxl_1.3.1           backports_1.2.1        Hmisc_4.4-2           
  [4] splines_4.0.3          crosstalk_1.1.1        BiocParallel_1.24.1   
  [7] digest_0.6.27          foreach_1.5.1          htmltools_0.5.1.1     
 [10] GO.db_3.12.1           checkmate_2.0.0        magrittr_2.0.1        
 [13] memoise_2.0.0          cluster_2.1.1          doParallel_1.0.16     
 [16] annotate_1.68.0        modelr_0.1.8           jpeg_0.1-8.1          
 [19] colorspace_2.0-0       blob_1.2.1             rvest_0.3.6           
 [22] haven_2.3.1            xfun_0.21              crayon_1.4.1          
 [25] RCurl_1.98-1.2         jsonlite_1.7.2         genefilter_1.72.1     
 [28] impute_1.64.0          survival_3.2-7         iterators_1.0.13      
 [31] glue_1.4.2             gtable_0.3.0           zlibbioc_1.36.0       
 [34] XVector_0.30.0         DelayedArray_0.16.3    scales_1.1.1          
 [37] DBI_1.1.1              Rcpp_1.0.6             htmlTable_2.1.0       
 [40] xtable_1.8-4           gridGraphics_0.5-1     foreign_0.8-81        
 [43] bit_4.0.4              preprocessCore_1.52.1  Formula_1.2-4         
 [46] htmlwidgets_1.5.3      httr_1.4.2             ellipsis_0.3.1        
 [49] farver_2.0.3           pkgconfig_2.0.3        XML_3.99-0.5          
 [52] nnet_7.3-15            dbplyr_2.1.0           locfit_1.5-9.4        
 [55] labeling_0.4.2         tidyselect_1.1.0       rlang_1.0.2           
 [58] later_1.1.0.1          munsell_0.5.0          cellranger_1.1.0      
 [61] tools_4.0.3            cachem_1.0.4           cli_2.3.0             
 [64] generics_0.1.0         RSQLite_2.2.3          broom_0.7.4           
 [67] evaluate_0.14          fastmap_1.1.0          yaml_2.2.1            
 [70] knitr_1.31             bit64_4.0.5            fs_1.5.0              
 [73] caTools_1.18.1         whisker_0.4            xml2_1.3.2            
 [76] compiler_4.0.3         rstudioapi_0.13        curl_4.3              
 [79] png_0.1-7              reprex_1.0.0           geneplotter_1.68.0    
 [82] stringi_1.5.3          highr_0.8              lattice_0.20-41       
 [85] Matrix_1.3-2           vctrs_0.3.6            pillar_1.4.7          
 [88] lifecycle_1.0.0        BiocManager_1.30.10    bitops_1.0-6          
 [91] httpuv_1.5.5           R6_2.5.0               latticeExtra_0.6-29   
 [94] promises_1.2.0.1       gridExtra_2.3          KernSmooth_2.23-18    
 [97] codetools_0.2-18       gtools_3.8.2           assertthat_0.2.1      
[100] rprojroot_2.0.2        rjson_0.2.20           withr_2.4.1           
[103] GenomeInfoDbData_1.2.4 hms_1.0.0              rpart_4.1-15          
[106] rmarkdown_2.6          rvcheck_0.1.8          git2r_0.28.0          
[109] lubridate_1.7.9.2      base64enc_0.1-3       

---

This R Markdown site was created with workflowr

---