...
We've done this several times before, so you should be able to come up with the full command lines if you refer back to the original lesson.
| Warning |
|---|
Be careful we are now mapping single-end reads, so you may have to look at the bowtie help to figure out how to do that! |
...
| Expand | ||||||||
|---|---|---|---|---|---|---|---|---|
| ||||||||
Now create a
Remember that there are 12 processors per node on Lonestar, so we choose to use 3 for each of the 4 jobs with the Create the launcher script and run it:
|
...
| Code Block |
|---|
head gene_counts.gff |
| Warning |
Optional: HTseq
HTseq is another tool to count reads. bedtools has many many useful functions, and counting reads is just one of them. In contrast, HTseq is a specialized utility for counting reads, and it does not have many functions other than that. HTseq is very slow and you need to run multiple command lines in order to do the same job as what bedtools multicov did. Why do we learn this? Well, you may want to care about reads mapped on intersection when you count reads. Please take a look at this page, and if this sophisticated counting method looks useful for you, use HTseq. Otherwise, use bedtools.
| Code Block |
|---|
grep "^NC_017544" NC_017544.1.gff > count_ref.gff
samtools view SRR034450.bam | htseq-count -m intersection-nonempty -t gene -i ID - count_ref.gff > count1.gff
samtools view SRR034451.bam | htseq-count -m intersection-nonempty -t gene -i ID - count_ref.gff > count2.gff
samtools view SRR034452.bam | htseq-count -m intersection-nonempty -t gene -i ID - count_ref.gff > count3.gff
samtools view SRR034453.bam | htseq-count -m intersection-nonempty -t gene -i ID - count_ref.gff > count4.gff
join count1.gff count2.gff | join - count3.gff | join - count4.gff > gene_counts_HTseq.gff
#if you have many samples, use for-loop and join
|
gene_counts_HTseq.gff has 5 more lines than gene_counts.gff. Check out the last 5 lines. They are basic statistics.
| Code Block |
|---|
wc -l gene_counts_HTseq.gff
tail gene_counts_HTseq.gff
|
The basic statistics (last 5 lines) is useful to know, but should be removed to use it as a input file for DEGseq
| Code Block |
|---|
head -2910 gene_counts_HTseq.gff > gene_counts_HTseq.tab
|
Finally, gene_counts_HTseq.tab is ready to use. HTseq-count is strand-specific in default. Therefore, read counts for each gene in gene_counts_HTseq.gff are approximately a half counts in gene_counts.gff for the corresponding gene.
Analyze differential gene expression
DESeq
Our data that is cluttered with a lot of extra columns and one column stuffed with tag=value information (including the gene names that we want!). Let's clean it up a bit before loading into R - which likes to work on simple tables. GFF are tab-delimited files.
We can do this cleanup many ways, but a quick one is to use the Unix string editor sed. This command replaces the entire beginning of the line up to locus_tag= with nothing (that is, it deletes it). This conveniently leaves us with just the locus_tag and the columns of read counts in each gene. If you were writing a real pipeline, you would probably want to use a Perl or Python script that would check to be sure that each line had the locus_tag (they do), among other things.
| Code Block | ||
|---|---|---|
| ||
head gene_counts.gff
sed 's/^.*locus_tag=//' gene_counts.gff > gene_counts.tab
|
After it has run, take a peek at the new file:
| Code Block |
|---|
head gene_counts.tab
|
| Warning | ||||
|---|---|---|---|---|
| ||||
Do not copy the
Is the same as:
It's ok to copy across the multiple lines and paste into R as long as you get all the way to the closing parenthesis. |
...
| Code Block | ||
|---|---|---|
| ||
login1$ R
...
> library("edgeR")
> counts = read.delim("gene_counts.tab", header=F, row.names=1)
> colnames(counts) = c("wt1", "mut1", "wt2", "mut2")
> head(counts)
> group <- factor(c("wt","mut","wt","mut"))
> dge = DGEList(counts=counts,group=group)
> dge <- estimateCommonDisp(dge)
> dge <- estimateTagwiseDisp(dge)
> et <- exactTest(dge)
> etp <- topTags(et, n=100000)
> etp$table$logFC = -etp$table$logFC
> pdf("edgeR-MA-plot.pdf")
> plot(
etp$table$logCPM,
etp$table$logFC,
xlim=c(-3, 20), ylim=c(-12, 12), pch=20, cex=.3,
col = ifelse( etp$table$FDR < .1, "red", "black" ) )
> dev.off()
> write.csv(etp$table, "edgeR-wt-vs-mut.csv")
> q()
Save workspace image? [y/n/c]: n
login1$ head edgeR-wt-vs-mut.csv
|
Note that the "FC" fold change calculated is initially the reverse of that for the DESeq example for the output here. It is wt relative to mut. To fix this, we put a negative in there for the log fold change.
...
- In an actual RNAseq analysis, you might want to trim stray adaptor sequences from your data using a tool like the FASTX-Toolkit, FAR, or cutadapt before aligningthe tools discussed in Evaluating your raw sequencing data.
- You can get a lot more information from RNAseq data than you could from a microarray experiment. You can map transcriptional start sites, areas of unexpected transcription, splice sites, etc. - all because you have full sequence information that we have barely used in this example.
- You can call variants from mapped RNAseq data, just be aware that many regions will have no coverage (because they are not expressed as RNA).
...