DESeq help

[emalandrak]

Hello everybody
I am new in RNAseq analysis, so I have some trouble trying to study differential gene expression. DEseq manual is really informative to work with but my question is about creating the input file. I found HT count but it needs a gff file and a SAM file with the aligned reads. If there is not a reference genome or transcriptome to download a gff file which is the way to construct it? I suppose by changing a tab delimited file from blastx but which exactly is the format for a gff as an input to HT count? Is there any automated software to annotate the constructed transcripts (output of e.g. Oases or Trinity) and create a gff?
Thanks in advance.

[simonandrews]

It might be easier to say what you do have and then work from there. Ultimately DESeq just needs a matrix of count values for different genes and it doesn't really care where they come from, so if you have count data already you can just make a matrix out of it and start from that point in DESeq.

[emalandrak]

DE seq

Dear Simon
Thank you very much for your reply.
The fact is that I don't have the count data. I have the BAM files from tophat, where the reads are mapped on the contructed transcriptome from Trinity. The problem is how to produce the count table for DESeq. I tryed HTseq but I had problems on installing it (and also needs a gff). So i need the way to produce the count table for DEseq

[simonandrews]

If the reads are mapped to assembled transcripts already then your count table is simply going to be the number of hits to each different accession in the set of transcripts. You just need to collate the number of hits to each different sequence id.

There may well be a program around which collates this already, but if not then it's a very small script to generate this data from a set of BAM files. If no one knows of a pre-built solution for this and you're not confident in having a go yourself then I can put something together which does this pretty quickly.

[simonandrews]

Having thought about this I found I did have a script which almost does what you want. This was something we used for miRNA mapping, but it's essentially the same idea you have in that you're mapping to a library of different sequences. The script below simply collates the hits to each different sequence in a set of SAM files (you could pipe the open of your bam files through "samtools view" to make it work from BAM files.

#!/usr/bin/perl
use warnings;
use strict;

my @sam_files = <*sam>;

my %counts;

for my $index (0..$#sam_files) {
    my $sam_file = $sam_files[$index];

    warn "Processing $sam_file\n";

    open (IN,$sam_file) or die $!;

    while (<IN>) {
	next if (/^@/);

	my (undef,undef,$mirna) = split(/\t/);
	next if ($mirna eq '*');

	$counts{$mirna}->[$index]++;
    }
}

open (OUT,'>','summarised_mirna_counts.txt') or die $!;

print OUT join("\t",("miRNA",@sam_files)),"\n";

foreach my $mirna (keys %counts) {

    my @counts = @{$counts{$mirna}};

    for (0..$#counts) {
	$counts[$_] = 0 unless ($counts[$_]);
    }

    for ($#counts+1..$#sam_files) {
	$counts[$_] = 0;
    }

    print OUT join("\t",($mirna,@counts)),"\n";
} 


close OUT or die $!;

[twotwo]

Hi,
I used HTseq from python. First you need to install python, then HTseq.... Just follow the instruction.

[emalandrak]

DEseq help

Dear Simon
Thank you very much for your script.
But I can't understand the script, since I have no background on programming. Anyway, I ll read some Pearl and I will give it a try. Thank you very much.

[emalandrak]

DEseq help

Thank you twotwo

I will give it a try again....