Hi,

I analyzing RNA-seq with two samples: control and treatment. I use the following sequence of software: sortmerna -> hisat2 -> DESeq2. This is my R code for DESeq2:

library( "GenomicFeatures" )
hse <- makeTxDbFromGFF( "~/Homo_sapiens.GRCh38.92.gtf", format="gtf" )
trans <- transcriptsBy( hse, by="gene" )

fls <- list.files( "./", pattern="bam$", full=TRUE )

library( "Rsamtools" )
bamLst <- BamFileList( fls, yieldSize=100000 )

library( "GenomicAlignments" )
se <- summarizeOverlaps( trans, bamLst,
mode="Union",
singleEnd=FALSE,
ignore.strand=TRUE,
fragments=TRUE )

library( "DESeq2" )

sampleInfo <- read.csv( "sample.csv" )
sampleInfo <- DataFrame( sampleInfo )
seIdx <- match(colnames(se), sampleInfo$run)

colData(se) <- cbind( colData(se), sampleInfo)

ddsFull <- DESeqDataSet( se, design= ~run )

countdata <- assay(se)
coldata <- colData (se)


ddsFullCountTable <- DESeqDataSetFromMatrix(
countData = countdata,
colData = coldata,
design = ~run)

dds <- ddsFullCountTable

dds <- DESeq(dds)
res <- results( dds )

write.csv( as.data.frame(res), file="results.csv" )

I wouldn’t recommend counting in transcripts unless you want intronic and incompatibly spliced reads to count for a gene. I’d recommend using the tximport pipeline, htseq-count or featureCounts with exon-by-gene features.

Otherwise, there isn’t a reason that baseMean would be rounded to 0.5 — try using makeExampleDESeqDataSet()