Estimating Gene Count Exercise - Revision history

Rf at 20:45, 11 May 2017

2017-05-11T20:45:33Z

Rf: Rf moved page I2rda egcx to Estimating Gene Count Exercise over redirect

2017-05-04T13:40:50Z

Rf moved page I2rda egcx to Estimating Gene Count Exercise over redirect

Rf: Rf moved page Estimating Gene Count Exercise to I2rda egcx

2017-05-04T12:54:32Z

Rf moved page Estimating Gene Count Exercise to I2rda egcx

Rf: s2

2017-04-08T22:32:09Z

Rf: s2

2017-04-08T22:31:24Z

Rf: s2

2017-04-08T22:30:12Z

Rf: s2

2017-04-08T22:26:08Z

Rf: Rf moved page Estimating Gen Count Exercise to Estimating Gene Count Exercise

2017-04-08T22:22:01Z

Rf moved page Estimating Gen Count Exercise to Estimating Gene Count Exercise

Rf: s1

2017-04-08T22:20:25Z

New page

= Aims =

In this part you will learn to:
generate count tables
You will use the following tools, which have been pre-installed on our bioinformatics training server at the University of
St Andrews:
HTSeq v0.6.1p1: http://www-huber.embl.de/users/anders/HTSeq/doc/count.html
The data set you'll be using is downloaded from ENA (http://www.ebi.ac.uk/ena/data/view/SRP019027). The reads belong to
samples SRR769314 and SRR769316. The data set is tailored with respect to the time allocated for the workshop. Reads were alig
ned to the first 20 Mb of chromosome 19 of the mouse reference genome (GRCm38/mm10) using TopHat and duplicates marked
using Picard MarkDuplicates.
You will use the following files:
SRR769314_duplicates_marked.bam: aligned reads
SRR769316_duplicates_marked.bam: aligned reads
mm10_chr19-1-20000000_Ensembl.gtf: Ensembl mouse gene models

Type text like this in the terminal at the $ command prompt, then press the
[Enter] key to run the command.

Data
The data is available in the directory 07_Estimating_gene_count:

cd /home/training/Data/07_Estimating_gene_count

Generate count tables
Generate a count table for sample SRR769314 using HTSeq htseq-count:

htseq-count -f bam -r pos -t exon \
-s no -m intersection-nonempty \
SRR769314_duplicates_marked.bam \
Reference/mm10_chr19-1-20000000_Ensembl.gtf \
> SRR769314_duplicates_marked_htseq.count

<ins>where</ins>:
* <code>-f</code>: format of the input data (default: sam).
* <code>-r</code>: for paired-end data, the alignment have to be sorted either by read name or by alignment position (default: name).
* <code>-t</code>: feature type (3rd column in GFF/GTF file) to be used, all features of other type are ignored (default: exon).
* <code>-s</code>: whether the data is from a strand-specific assay (default: yes).
* <code>-m</code>: mode to handle reads overlapping more than one feature (default: union).

Repeat for sample SRR769316:

htseq-count -f bam -r pos -t exon \
-s no -m intersection-nonempty \
SRR769316_duplicates_marked.bam \
Reference/mm10_chr19-1-20000000_Ensembl.gtf \
> SRR769316_duplicates_marked_htseq.count

Have a look at the start and end of the results files:

head SRR769314_duplicates_marked_htseq.count
tail SRR769314_duplicates_marked_htseq.count
head SRR769316_duplicates_marked_htseq.count
tail SRR769316_duplicates_marked_htseq.count

Count the proportion of reads that fall in annotated genes, using a small AWK script:

awk '{if(/^__/){excluded+=$2}else{included+=$2}} \
END{printf "%s,%s,%.2f%%\n", included, excluded,\
included/(included+excluded)*100}' \
SRR769314_duplicates_marked_htseq.count
awk '{if(/^__/){excluded+=$2}else{included+=$2}} \
END{printf "%s,%s,%.2f%%\n", included, excluded, \
included/(included+excluded)*100}' \
SRR769316_duplicates_marked_htseq.count

which means:
for each line in the count file:
if the line contains the pattern __ at the start of the line
increase the value of the variable excluded with the value in the second field of the line
else increase the value of the variable included with the value in the second field of the line
then:
print the value of included, excluded and included/(included+excluded)*100, as a string, a string and a floating point number with two decimal place, respectively, followed by a line break

* What proportion of read are excluded from further analysis?
* What is the main reason for reads being excluded?

Run HTSeq htseq-count again changing the mode from intersection-nonempty to union.
* What difference does this make?

@@ Line 16: / Line 16: @@
 = Data =
-The data is available in the directory <code>07_Estimating_gene_count</code>:
+The data is available in the directory <code>05_Estimating_gene_count</code>:
- cd /home/training/Data/07_Estimating_gene_count
+  cd ~/i2rda_data/05_Estimating_gene_count
 = Generate count tables =
 Generate a count table for sample SRR769314 using HTSeq's <code>htseq-count</code>:
-  htseq-count -f bam -r pos -t exon \
+  htseq-count -f bam -r pos -t exon -s no -m intersection-nonempty SRR769314_duplicates_marked.bam Reference_files/mm10_chr19-1-20000000_Ensembl.gtf > SRR769314_duplicates_marked_htseq.count
 <ins>where</ins>:
@@ Line 36: / Line 32: @@
 * <code>-m</code>: mode to handle reads overlapping more than one feature (default: union).
-Repeat for sample SRR769316:
+We want to repeat this command for sample <code>SRR769316</code>. Use the history command and the command-line shortcuts to do it.
- htseq-count -f bam -r pos -t exon \
+<ins>Question</ins>:
- -s no -m intersection-nonempty \
+* Can you imagine a situation where this could cause problems?
- SRR769316_duplicates_marked.bam \
+* How might you remedy them?
-Have a look at the start and end of the results files:
+Have a look at the start and end of the results files, we'll use the replacement short cut <code>^this^that</code>.
   head SRR769314_duplicates_marked_htseq.count
   tail SRR769314_duplicates_marked_htseq.count
-  head SRR769316_duplicates_marked_htseq.count
+  ^14^16
-Count the proportion of reads that fall in annotated genes, using a small AWK script:
+<ins>Question</ins>:
-  awk '{if(/^__/){excluded+=$2}else{included+=$2}} \
+* Two columns, first is gene name, then a count value.
- END{printf "%s,%s,%.2f%%\n", included, excluded,\
+* some metadata at the end of the file, indicated by leading underscores.
- included/(included+excluded)*100}' \
- SRR769314_duplicates_marked_htseq.count
+We're about to embark on an <code>awk</code> writing exercise, let's start with an easy one. We want to put the metadata in perspective. We want to know the total number of reads:
- awk '{if(/^__/){excluded+=$2}else{included+=$2}} \
+<ins>How this AWK command works</ins>:
   for each line in the count file:
   if the line contains the pattern __ at the start of the line
-  increase the value of the variable excluded with the value in the second field of the line
+  increase the value of the variable OUT with the value in the second field of the line
-  else increase the value of the variable included with the value in the second field of the line
+  else increase the value of the variable IN with the value in the second field of the line
   then:
-  print the value of included, excluded and included/(included+excluded)*100, as a string, a string and a floating point number with two decimal place, respectively, followed by a line break
+  print the value of IN, OUT and IN/(IN+OUT)*100, as a string, a string and a floating point number with two decimal place, respectively, followed by a line break
 * What proportion of reads are excluded from further analysis?
 * What is the main reason for reads being excluded?

← Older revision	Revision as of 13:40, 4 May 2017
(No difference)

← Older revision	Revision as of 12:54, 4 May 2017
(No difference)

@@ Line 71: / Line 71: @@
   print the value of included, excluded and included/(included+excluded)*100, as a string, a string and a floating point number with two decimal place, respectively, followed by a line break
-* What proportion of read are excluded from further analysis?
+* What proportion of reads are excluded from further analysis?
 * What is the main reason for reads being excluded?
 Run HTSeq htseq-count again changing the mode from intersection-nonempty to union.
 * What difference does this make?

@@ Line 4: / Line 4: @@
 * generate count tables
-HTSeq v0.6.1p1: http://www-huber.embl.de/users/anders/HTSeq/doc/count.html
+The tool you will use is:
 The data set you'll be using is downloaded from ENA (http://www.ebi.ac.uk/ena/data/view/SRP019027). The reads belong to
-samples SRR769314 and SRR769316. The data set is tailored with respect to the time allocated for the workshop. Reads were alig
+samples <code>SRR769314</code> and <code>SRR769316</code>. The data set is tailored with respect to the time allocated for the exercise. Reads were aligned to the first 20 Mb of chromosome 19 of the mouse reference genome (<code>GRCm38/mm10</code>) using <code>TopHat</code> and duplicates marked using Picard- tools' <code>MarkDuplicates</code>.
-ned to the first 20 Mb of chromosome 19 of the mouse reference genome (GRCm38/mm10) using TopHat and duplicates marked
 You will use the following files:
-SRR769314_duplicates_marked.bam: aligned reads
+* <code>SRR769314_duplicates_marked.bam</code>: aligned reads
-SRR769316_duplicates_marked.bam: aligned reads
+* <code>SRR769316_duplicates_marked.bam</code>: aligned reads
-mm10_chr19-1-20000000_Ensembl.gtf: Ensembl mouse gene models
+* <code>mm10_chr19-1-20000000_Ensembl.gtf</code>: Ensembl mouse gene models
 = Data =

← Older revision	Revision as of 22:22, 8 April 2017
(No difference)