DGE of RNAseq data using DeSEQ2 - need help with input file

Question

[Wilhelm]

Hey all,

I'm new to RNAseq and am having a little trouble getting my data in a format usable for DGE programs. I tried using cufflinks with my sorted indexed bam files and a .gff file but it keeps throwing back a "segmentation fault (core dumped)" error after loading reference annotation; I think this is a hardware error so I have no idea how to get around it.

My workflow thus far:

• Illumina miSEQ 75bp PE reads of untreated vs treated (2 biological replicates each condition); working on bacteria
• assembled and aligned using bowtie2 against my reference genome
• converted to bam with samtools; sorted and indexed
• run through bedtools multicov for 1 output .bed file with counts for each gene and filtered for >Q30

So I end up with a .bed file with my 4 counts (untreated1 untreated 2 treated 1 treated 2).

I want to use deSEQ2 and edgeR, but have no experience with using R at all.

Would anyone be able to give some advice on a basic script workflow to use for deSEQ2 and edgeR? I've been perusing forums for a day or so now and reading to no avail.

Cheers

[Martombo]

the easiest thing you can do is probably to use htseq-count instead of bedtools to get the gene counts. then you will find a function in deseq2 to directly import the results. you can find a description on deseq2 guide.

[Wilhelm]

I should mention I've gotten as far as this:

directory<-"/deSEQ2"
sampleFiles <- grep("",list.files(directory),value=TRUE)
sampleFiles
[1] "GSH1.bed" "GSH2.bed" "WT1.bed"  "WT2.bed" 
sampleCondition<-c("2mM GSH","2mM GSH","WT","WT")
sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, condition=sampleCondition)
sampleTable
  sampleName fileName condition
1   GSH1.bed GSH1.bed   2mM GSH
2   GSH2.bed GSH2.bed   2mM GSH
3    WT1.bed  WT1.bed        WT
4    WT2.bed  WT2.bed        WT

I have a count matrix (all 4 counts in a bed file), and also four individual bed files; my question is how do I get my counts to load into my design?

[Martombo]

if you really want to use bedtools output instead, you should read in R your files with a function like read.table(). then you can merge the different tables and transform the counts to a matrix, which can be used as input for deseq2. but it's probably both safer and easier to use htseq-count.