Skip to content

Commit

Permalink
Merge pull request #863 from maxplanck-ie/develop
Browse files Browse the repository at this point in the history
Develop
  • Loading branch information
katsikora authored Dec 13, 2022
2 parents 46e5465 + 4cca146 commit 2b3fa05
Show file tree
Hide file tree
Showing 22 changed files with 546 additions and 110 deletions.
43 changes: 23 additions & 20 deletions .ci_stuff/test_dag.sh
Original file line number Diff line number Diff line change
Expand Up @@ -184,47 +184,50 @@ if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 333 ]; then exit 1 ; fi

# mRNA-seq
WC=`mRNA-seq -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 882 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 879 ]; then exit 1 ; fi
WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 892 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 889 ]; then exit 1 ; fi
WC=`mRNA-seq -i PE_input -o output --rMats --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 908 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 905 ]; then exit 1 ; fi
WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 602 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 599 ]; then exit 1 ; fi
WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment,deepTools_qc" --trim .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 948 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 945 ]; then exit 1 ; fi
WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment-free,deepTools_qc" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1005 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1006 ]; then exit 1 ; fi
WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment,deepTools_qc" --bcExtract --trim .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 916 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 913 ]; then exit 1 ; fi
WC=`mRNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment,deepTools_qc" --bcExtract --UMIDedup --trim .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 966 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 963 ]; then exit 1 ; fi
WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 807 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 804 ]; then exit 1 ; fi
WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 526 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 523 ]; then exit 1 ; fi
WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment,deepTools_qc" --trim .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 863 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 860 ]; then exit 1 ; fi
WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" -m "alignment-free,deepTools_qc" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 920 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 921 ]; then exit 1 ; fi
WC=`mRNA-seq -i SE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --trim --fastqc .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 975 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 972 ]; then exit 1 ; fi
WC=`mRNA-seq -i BAM_input/filtered_bam -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --fromBAM .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 659 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 657 ]; then exit 1 ; fi
#multiple comparison groups
WC=`mRNA-seq --mode alignment,alignment-free -i PE_input -o output --rMats --sampleSheet .ci_stuff/test_sampleSheet_multiComp.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 867 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 869 ]; then exit 1 ; fi
#allelic
WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1,strain2 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1360 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1357 ]; then exit 1 ; fi
WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i allelic_BAM_input/filtered_bam --fromBAM --bamExt '.filtered.bam' -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --SNPfile allelic_input/snpfile.txt --NMaskedIndex allelic_input/Ngenome .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1118 ]; then exit 1 ; fi
WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1,strain2 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1370 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1367 ]; then exit 1 ; fi
WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --SNPfile allelic_input/snpfile.txt --NMaskedIndex allelic_input/Ngenome .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1353 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1350 ]; then exit 1 ; fi
WC=`mRNA-seq -m allelic-mapping,deepTools_qc -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1370 ]; then exit 1 ; fi
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1367 ]; then exit 1 ; fi
WC=`mRNA-seq -m allelic-mapping,deepTools_qc,alignment-free -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" --VCFfile allelic_input/file.vcf.gz --strains strain1 .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 1803 ]; then exit 1 ; fi

# noncoding-RNA-seq
WC=`noncoding-RNA-seq -i PE_input -o output --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
if [ ${PIPESTATUS[0]} -ne 0 ] || [ $WC -ne 717 ]; then exit 1 ; fi
WC=`noncoding-RNA-seq -i PE_input -o output --sampleSheet .ci_stuff/test_sampleSheet.tsv --snakemakeOptions " --dryrun --conda-prefix /tmp" .ci_stuff/organism.yaml | tee >(cat 1>&2) | grep -v "Conda environment" | sed '/^\s*$/d' | wc -l`
Expand Down
2 changes: 1 addition & 1 deletion conda-recipe/meta.yaml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
package:
name: snakepipes
version: 2.6.1
version: 2.7.0

source:
path: ../
Expand Down
9 changes: 9 additions & 0 deletions docs/content/News.rst
100755 → 100644
Original file line number Diff line number Diff line change
@@ -1,11 +1,19 @@
snakePipes News
===============


snakePipes 2.7.0
----------------

* Added the allelic version of Salmon-based transcript quantitation to mRNA-seq workflow. Will be run if *both* 'allelic-mapping' and 'alignment-free' modes are specified. An allelic version of sleuth will be run, if sample sheet is provided, as well as DESeq2 on allelic Salmon counts.


snakePipes 2.6.1
----------------

* Capped tabulate version as 0.9.0 breaks snakemake


snakePipes 2.6.0
----------------

Expand All @@ -19,6 +27,7 @@ snakePipes 2.6.0
* Fixed a couple of issues in the ATAC-seq workflow after sofware versions update.
* Fixed genome size conversion to string.


snakePipes 2.5.4
----------------

Expand Down
2 changes: 1 addition & 1 deletion snakePipes/__init__.py
Original file line number Diff line number Diff line change
@@ -1 +1 @@
__version__ = '2.6.1'
__version__ = '2.7.0'
1 change: 1 addition & 0 deletions snakePipes/common_functions.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ def set_env_yamls():
return {'CONDA_SHARED_ENV': 'envs/shared.yaml',
'CONDA_CREATE_INDEX_ENV': 'envs/createIndices.yaml',
'CONDA_RNASEQ_ENV': 'envs/rna_seq.yaml',
'CONDA_SLEUTH_ENV': 'envs/sleuth.yaml',
'CONDA_RMATS_ENV': 'envs/rMats.yaml',
'CONDA_scRNASEQ_ENV': 'envs/sc_rna_seq.yaml',
'CONDA_seurat_ENV': 'envs/sc_rna_seq_seurat.yaml',
Expand Down
24 changes: 15 additions & 9 deletions snakePipes/shared/rscripts/DESeq2.R
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ if(isTRUE(tximport)) {
tx2gene <- read.delim(tx2gene_file, header = FALSE)
tx2gene <- tx2gene[c(1,2)]
# check setup table and import
countdata <- checktable(sampleSheet = sampleInfo, salmon_dir = dirname(countFilePath), tx2gene_annot = tx2gene)
countdata <- checktable(sampleSheet = sampleInfo, salmon_dir = dirname(countFilePath), tx2gene_annot = tx2gene, alleleSpecific = allelic_info)
} else {
sampleInfo$name <- make.names(sampleInfo$name)
rownames(sampleInfo)<-sampleInfo$name
Expand All @@ -93,19 +93,23 @@ if(isTRUE(tximport)) {
## ~~~~~~~ 3. run DESeq wrapper ~~~~~~~~
#in case of the allelic-specific workflow, allow for 1 condition and skip deseq2 basic in this case
if(length(unique(sampleInfo$condition))>1){
seqout <- DESeq_basic(countdata, coldata = sampleInfo, fdr = fdr, alleleSpecific = allelic_info, from_salmon = tximport)
if(tximport & allelic_info){
message("Detected allelic Salmon counts. Skipping DESeq_basic.")
}else{
seqout <- DESeq_basic(countdata, coldata = sampleInfo, fdr = fdr, alleleSpecific = allelic_info, from_salmon = tximport)

DESeq_writeOutput(DEseqout = seqout,
DESeq_writeOutput(DEseqout = seqout,
fdr = fdr, outprefix = "DEseq_basic",
geneNamesFile = geneNamesFilePath)
}
}
#DESeq_downstream(DEseqout = seqout, countdata, sampleInfo,
# fdr = fdr, outprefix = "DEseq_basic", heatmap_topN = topN,
# geneNamesFile = geneNamesFilePath)

## Run allele-sepecific DESeq wrapper (if asked for)
if (isTRUE(allelic_info)) {
seqout_allelic <- DESeq_allelic(countdata, coldata = sampleInfo, fdr = fdr)
seqout_allelic <- DESeq_allelic(countdata, coldata = sampleInfo, fdr = fdr, from_salmon=tximport)

DESeq_writeOutput(DEseqout = seqout_allelic,
fdr = fdr, outprefix = "DEseq_allelic",
Expand All @@ -132,19 +136,21 @@ write.bibtex(bib, file = 'citations.bib')
file.copy(rmdTemplate, to = 'DESeq2_report_basic.Rmd')

if(length(unique(sampleInfo$condition))>1){
outprefix = "DEseq_basic"
cite_options(citation_format = "text",style = "html",cite.style = "numeric",hyperlink = TRUE)
render('DESeq2_report_basic.Rmd',
if(!tximport & !allelic_info){
outprefix = "DEseq_basic"
cite_options(citation_format = "text",style = "html",cite.style = "numeric",hyperlink = TRUE)
render('DESeq2_report_basic.Rmd',
output_format = "html_document",
clean = TRUE,
params = list(
DEseqoutRdata = paste0(outprefix, "_DESeq.Rdata"),
ddr.df = paste0(outprefix, "_DEresults.tsv"),
countdata = countFilePath,
coldata = sampleInfo,
fdr = 0.05,
fdr = fdr,
heatmap_topN = 20,
geneNamesFile = geneNamesFilePath))
}
}

if (isTRUE(allelic_info)) {
Expand All @@ -158,7 +164,7 @@ if (isTRUE(allelic_info)) {
ddr.df = paste0(outprefix, "_DEresults.tsv"),
countdata = countFilePath,
coldata = sampleInfo,
fdr = 0.05,
fdr = fdr,
heatmap_topN = 20,
geneNamesFile = geneNamesFilePath))
}
Expand Down
2 changes: 1 addition & 1 deletion snakePipes/shared/rscripts/DESeq2Report.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ load(params$DEseqoutRdata) ## dds and ddr
ddr.df <- params$ddr.df #_DEresults.tsv file
countdata <- params$countdata
coldata <- params$coldata
fdr <- params$fdr
fdr <- as.numeric(params$fdr)
heatmap_topN <- params$heatmap_topN
geneNamesFile <- params$geneNamesFile
Expand Down
73 changes: 48 additions & 25 deletions snakePipes/shared/rscripts/DE_functions.R
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,13 @@ checktable <- function(countdata = NA, sampleSheet = NA, alleleSpecific = FALSE,
## check files
if(!is.na(salmon_dir)) {

# mode = Salmon : check whether salmon output files exist in Salmon dir
files <- paste0(salmon_dir,"/",sampleSheet$name,".quant.sf")
#mode = Salmon : check whether salmon output files exist in Salmon dir
if(alleleSpecific){
files <- c(paste0(salmon_dir,"/",sampleSheet$name,".genome1.quant.sf"),paste0(salmon_dir,"/",sampleSheet$name,".genome2.quant.sf"))
}else{
files <- paste0(salmon_dir,"/",sampleSheet$name,".quant.sf")
}
#files<-dir(salmon_dir,pattern=".quant.sf",full.names=TRUE)
filecheck <- file.exists(files)
if(!(all(filecheck == TRUE))) {
cat("Error! The following File(s) don't exist : ")
Expand Down Expand Up @@ -63,7 +68,7 @@ checktable <- function(countdata = NA, sampleSheet = NA, alleleSpecific = FALSE,
}
}
}
if(all(is.integer(countdata))){
if(all(is.integer(countdata)) | !is.na(salmon_dir)){
print("All countdata is integer.")
}else{
print("Non-integer counts detected. The data will be rounded, as this is well within the expected sampling variation of a technical replicate.")
Expand All @@ -88,25 +93,26 @@ DESeq_basic <- function(countdata, coldata, fdr, alleleSpecific = FALSE, from_sa
cnames.sub<-unique(colnames(coldata)[2:which(colnames(coldata) %in% "condition")])
d<-as.formula(noquote(paste0("~",paste(cnames.sub,collapse="+"))))


# Normal DESeq
print("Performing basic DESeq: test vs control")
if(isTRUE(from_salmon)) {
print("Using input from tximport")
dds <- DESeq2::DESeqDataSetFromTximport(countdata,
colData = coldata, design =d)

} else {
print("Using input from count table")
if(isTRUE(alleleSpecific)) {
rnasamp <- dplyr::select(countdata, dplyr::ends_with("_all"))
rownames(coldata)<-colnames(rnasamp)
dds <- DESeq2::DESeqDataSetFromMatrix(countData = rnasamp,
colData = coldata, design =d)
} else {
dds <- DESeq2::DESeqDataSetFromMatrix(countData = countdata,
if(isTRUE(alleleSpecific)) {
rnasamp <- dplyr::select(countdata, dplyr::ends_with("_all"))
rownames(coldata)<-colnames(rnasamp)
countdata<-rnasamp
}

dds <- DESeq2::DESeqDataSetFromMatrix(countData = countdata,
colData = coldata, design =d)

}
}

if(length(size_factors) > 1) {
print("applying size factors")
print(size_factors)
Expand Down Expand Up @@ -135,33 +141,50 @@ DESeq_basic <- function(countdata, coldata, fdr, alleleSpecific = FALSE, from_sa
#'
#'

DESeq_allelic <- function(countdata, coldata, fdr) {
DESeq_allelic <- function(countdata, coldata, fdr, from_salmon=FALSE) {

# AlleleSpecific DEseq
print("Performing Allele-specific DESeq using Interaction design : Genome2 vs Genome1")
if(isTRUE(from_salmon)) {
# create alleleSpecific design matrix
coldata_allelic <- data.frame(name = colnames(as.data.frame(countdata$counts)),
allele = rep(c("genome1", "genome2"), nrow(coldata)),
condition = rep(coldata$condition, each = 2) )
rownames(coldata_allelic)<-colnames(as.data.frame(countdata$counts))
coldata_allelic$allele<-factor(coldata_allelic$allele,levels=c("genome1","genome2"))
coldata_allelic$condition<-factor(coldata_allelic$condition,levels=unique(coldata_allelic$condition))
print("Using input from tximport")
dds <- DESeq2::DESeqDataSetFromTximport(countdata,
colData = coldata_allelic, design =~1)

} else {
print("Using input from count table")
rnasamp <- dplyr::select(countdata, -dplyr::ends_with("_all"))

# create alleleSpecific design matrix
design <- data.frame(name = colnames(rnasamp),
coldata_allelic <- data.frame(name = colnames(rnasamp),
allele = rep(c("genome1", "genome2"), nrow(coldata)),
condition = rep(coldata$condition, each = 2) )
rownames(design)<-colnames(rnasamp)
rownames(coldata_allelic)<-colnames(rnasamp)
coldata_allelic$allele<-factor(coldata_allelic$allele,levels=c("genome1","genome2"))
coldata_allelic$condition<-factor(coldata_allelic$condition,levels=unique(coldata_allelic$condition))
dds <- DESeq2::DESeqDataSetFromMatrix(rnasamp, colData = coldata_allelic,
design = ~1)
rownames(dds) <- rownames(rnasamp)

}

# Run DESeq
if(length(unique(design$condition))>1){
dds <- DESeq2::DESeqDataSetFromMatrix(rnasamp, colData = design,
design = ~allele + condition + allele:condition)
rownames(dds) <- rownames(rnasamp)
if(length(unique(coldata_allelic$condition))>1){
DESeq2::design(dds) <- formula(~allele + condition + allele:condition)
dds <- DESeq2::DESeq(dds,betaPrior = FALSE)
ddr <- DESeq2::results(dds, name=paste0("allelegenome2.condition",unique(coldata$condition)[2]))
ddr_shrunk <- DESeq2::lfcShrink(dds,coef=paste0("allelegenome2.condition",unique(coldata$condition)[2]),type="apeglm",res=ddr)
} else {
dds <- DESeq2::DESeqDataSetFromMatrix(rnasamp, colData = design,
design = ~allele)
rownames(dds) <- rownames(rnasamp)
dds <- DESeq2::DESeq(dds,betaPrior = FALSE)
ddr <- DESeq2::results(dds, name="allele_genome2_vs_genome1")
ddr_shrunk <- DESeq2::lfcShrink(dds,coef="allele_genome2_vs_genome1",type="apeglm",res=ddr)
DESeq2::design(dds) <- formula(~allele)
dds <- DESeq2::DESeq(dds,betaPrior = FALSE)
ddr <- DESeq2::results(dds, name="allele_genome2_vs_genome1")
ddr_shrunk <- DESeq2::lfcShrink(dds,coef="allele_genome2_vs_genome1",type="apeglm",res=ddr)
}
output <- list(dds = dds, ddr = ddr, ddr_shrunk=ddr_shrunk)
return(output)
Expand Down
Loading

0 comments on commit 2b3fa05

Please sign in to comment.