Answers

To run the commands alongside the answers you must first have run all of the tutorial commands for that section and generated the output files they reference.

Tutorial sections

• Introducing the tutorial dataset
• Mapping RNA-Seq reads to the genome using HISAT2
• Visualising transcriptomes with IGV
• Transcript quantification with Kallisto
• Identifying differentially expressed genes with Sleuth
• Interpreting the results

Additional information

• Normalisation

Let's go to our data directory.

Introducing the tutorial dataset

Q1: Why is there more than one FASTQ file per sample?

There are 2 FASTQ files for each sample e.g. MT1_1.fastq and MT1_2.fastq. This is because this was paired-end sequence data.

With Illumina paired-end sequencing, you will have a fragment (hopefully longer than your reads) which is sequenced at both ends. This means that there will be a "left" read (sometimes known as r1) and its corresonding "right" read (sometimes known as r2). These are indicated by the /1 and /2 at the end of the read name in the _1.fastq and _2.fastq files respectively.

So, the left read header _@HS1808296:8:1101:1352:48181#4 in MT1_1.fastq has the /1 suffix:

@HS18_08296:8:1101:1352:48181#4/1 

And, the corresponding right read header in MT1_2.fastq has the /2 suffix:

@HS18_08296:8:1101:1352:48181#4/2

Q2: How many reads are there for MT1?

There are 2.5 million reads for MT1.

Ideally, you should count the reads in both files. This is because sometimes we have singletons (reads without a mate) after preprocessing steps such as trimming.

Our reads look like:

@HS18_08296:8:1101:1352:48181#4/2
ATCCGCCNANTTTNNNNATATAATTANNNGNAANNAANNNNAATNACANNNATTNNNNTAGNANNNGNNAGTNNACAAGGNTNNNNNNNAAAGNNNNNNN
+
9AABE8A!D!DFA!!!!EE@CFCD@B!!!D!F6!!EE!!!!EE4!F5E!!!BEA!!!!B@6!A!!!D!!FD'!!C+D@@*!B!!!!!!!B>DE!!!!!!!

With the FASTQ format there are four lines per read. So, as long as our files are not truncated we can count the number of lines and divide them by four.

So, we can then divide this by four to give us 1.25 million. We will get the same for our r2 reads:

So, we have 1.25 million read pairs or 2.5 million (1.25 x 2) reads for our MT1 sample.

You may also have thought "aha, I can use the @ which is at the start of the header"...

But, wait, there are 1,343,714 left reads and 1,250,000 right reads...can that be right...no.

Take a closer look at the quality scores on the fourth line...they also contain @. This is because the quality scores are Phred+33 encoded. For more information on quality score encoding see our Data Formats and QC tutorial or go here.

Instead, we can use the earlier information about the left and right read suffixes: /1 and /2. This can be with two commands:

Or, with only one command:

Note: Don't forget to use the -c option for grep to count the occurences and $ to make sure you're only looking for \1 or \2 at the end of the line.

Mapping RNA-Seq reads to the genome using HISAT2

Map, convert (SAM to BAM), sort and index using the reads from the MT2 sample.

You can do this with several, individual steps:

Or, you can do it in one step:

Or, you could write a for loop to run the command above for all five of your samples:

For more information on how this loop works, have a look at our Unix tutorial and Running commands on multiple samples.

Q1: How many index files were generated when you ran hisat2-build?

There are 8 HISAT2 index files for our reference genome.

Q2: What was the overall alignment rate for each of the MT samples (MT1 and MT2) to the reference genome?

The overall alignment rate for MT1 was 94.15% and 91.69% for MT2.

1250000 reads; of these:
  1250000 (100.00%) were paired; of these:
    105798 (8.46%) aligned concordantly 0 times
    329468 (26.36%) aligned concordantly exactly 1 time
    814734 (65.18%) aligned concordantly >1 times
    ----
    105798 pairs aligned concordantly 0 times; of these:
      1797 (1.70%) aligned discordantly 1 time
    ----
    104001 pairs aligned 0 times concordantly or discordantly; of these:
      208002 mates make up the pairs; of these:
        146250 (70.31%) aligned 0 times
        19845 (9.54%) aligned exactly 1 time
        41907 (20.15%) aligned >1 times
94.15% overall alignment rate


1250000 reads; of these:
  1250000 (100.00%) were paired; of these:
    139557 (11.16%) aligned concordantly 0 times
    483583 (38.69%) aligned concordantly exactly 1 time
    626860 (50.15%) aligned concordantly >1 times
    ----
    139557 pairs aligned concordantly 0 times; of these:
      4965 (3.56%) aligned discordantly 1 time
    ----
    134592 pairs aligned 0 times concordantly or discordantly; of these:
      269184 mates make up the pairs; of these:
        207729 (77.17%) aligned 0 times
        28836 (10.71%) aligned exactly 1 time
        32619 (12.12%) aligned >1 times
91.69% overall alignment rate

Note: If a read pair is concordantly aligned it means both reads in the pair align with the same chromosome/scaffold/contig, the reads are aligned in a proper orientation (typically ----> <----) and that the reads have an appropriate insert size.

Q3: How many MT1 and MT2 reads were not aligned to the reference genome?

146,250 reads (5.85%) and 207,729 reads (8.31%) did not align to the reference genome for MT1 and MT2 respectively.

Here is a brief summary of what the HISAT2 summary tells us for our MT2 sample and how we can tell which of the summary lines gives us this information:

• We have 1,250,000 read pairs or 2,500,000 reads (2 x 1,250,000 pairs)

   1250000 reads; of these:

• All of our reads (100%) are paired - i.e. no reads without their mate

  ``` 1250000 (100.00%) were paired; of these:

```

• 1,110,443 pairs (88.84%) align concordantly one (38.69%) or more (50.15%) times

   483583 (38.69%) aligned concordantly exactly 1 time
   626860 (50.15%) aligned concordantly >1 times

• 139,557 pairs (11.16%) or 279,114 reads (2 x 139,557) did not align concordantly anywhere in the genome

   139557 (11.16%) aligned concordantly 0 times

• Of those 139,557 pairs, 4,965 pairs align discordantly (3.56% of the 139,557 pairs)

   139557 pairs aligned concordantly 0 times; of these:
    4965 (3.56%) aligned discordantly 1 time

• This leave us with 134,592 pairs (139,557 - 4,965) where both reads in the pair do not align to the genome (concordantly or discordantly)

   134592 pairs aligned 0 times concordantly or discordantly; of these:

• Of those 269,184 reads (2 x 134,592) we have 61,455 reads (22.83%) which align to the genome without their mate

   269184 mates make up the pairs; of these:
   ...
   28836 (10.71%) aligned exactly 1 time
   32619 (12.12%) aligned >1 times

• Leaving us with a 91.69% overall alignment rate

  91.69% overall alignment rate

• That means 207,729 (8.31%) of the 2,500,000 reads (or 77.17% of the unaligned pairs) do not align anywhere in the genome

  207729 (77.17%) aligned 0 times

Visualising transcriptomes with IGV

Q1: How many CDS features are there in "PCHAS_1402500"?

There are 8 CDS features in PCHAS_1402500. You can get this in several ways:

Count the number of exons/CDS features in the gene annotation.

Count the number of CDS features in the GFF file.

First, get all of the CDS features for PCHAS_1402500.

Then count the number of the CDS features for PCHAS_1402500.

Q2: Does the RNA-seq mapping agree with the gene model in blue?

Yes. The peaks of the coverage tracks correspond to the annotated exon/CDS features.

Q3: Do you think this gene is differentially expressed and is looking at the coverage plots alone a reliable way to assess differential expression?

Possibly. But, you can't tell differential expression by the counts alone as there may be differences in the sequencing depths of the samples.

Transcript quantification with Kallisto

Use kallisto quant four more times, for the MT2 sample and the three SBP samples.

You can run the individual kallisto quant commands as you did for MT1 for each of the remaining samples:

Or, you can write a for loop which will run kallisto quant on all of the samples:

For more information on how this loop works, have a look at our Unix tutorial and Running commands on multiple samples.

Q1: What k-mer length was used to build the Kallisto index?

A k-mer length of 31 was used.

Look at the output from kallisto index:

[build] k-mer length: 31

Or, look for the -k or --kmer-size option in the kallisto index usage:

Q2: How many transcript sequences are there in PccAS_v3_transcripts.fa?

There are 5177 transcript sequences.

Look at the output from kallisto quant:

[index] number of targets: 5,177

Or, look for n_targets in one of the run_info.json files:

Or, you can run a grep on the transcript FASTA file and count the number of header lines:

Q3: What is the transcripts per million (TPM) value for PCHAS_1402500 in each of the samples?

You can look at each of the individual abundance files:

Or you can use a recursive grep:

Or you can use a loop:

For more information on how this loop works, have a look at our Unix tutorial and Running commands on multiple samples.

Q4: Do you think PCHAS_1402500 is differentially expressed?

Probably not. We would need to run statistical tests to really be sure though.

Identifying differentially expressed genes with Sleuth

Q1: Is our gene from earlier, PCHAS_1402500, significantly differentially expressed?

No.

Look at the transcript table.

And the transcript view.

Although this gene looked like it was differentially expressed from the plots in IGV, our test did not show it to be so (q-value > 0.05). This might be because some samples tended to have more reads, so based on raw read counts, genes generally look up-regulated in the SBP samples.

Alternatively, the reliability of only two biological replicates and the strength of the difference between the conditions was not sufficient to be statistically convincing. In the second case, increasing the number of biological replicates would give us more confidence about whether there really was a difference.

In this case, it was the lower number of reads mapping to MT samples that mislead us in the IGV view. Luckily, careful normalisation and appropriate use of statistics saved the day!

Interpreting the results

Q1: How many genes are more highly expressed in the SBP samples?

127. We can use awk to filter our Kallisto/sleuth results and wc -l to count the number of lines returned.

Q2: How many genes are more highly expressed in the MT samples?

169. We can use awk to filter our Kallisto/sleuth results and wc -l to count the number of lines returned.

Q3: Do you notice any particular genes that come up in the analysis?

Genes from the cir family are upregulated in the MT samples.

To get this we must first find out which genes (or gene descriptions) are seen most often in the genes which are more highly expressed in our SBP samples.

Then, we summarise the genes more highly expressed in our MT samples.

Perhaps the CIR proteins are interesting. There are only 2 CIR proteins upregulated in the SBP samples and 25 CIR in the MT samples.

The cir family is a large, malaria-specific gene family which had previously been proposed to be involved in immune evasion (Lawton et al., 2012). Here, however, we see many of these genes upregulated in a form of the parasite which seems to cause the immune system to better control the parasite. This suggests that these genes interact with the immune system in a subtler way, preventing the immune system from damaging the host.

Normalisation

How we got the information to help the questions

To answer the questions you needed the following for each sample:

• Number of reads assigned to PCHAS_1402500
• Length of exons in PCHAS_1402500 (bp)
• Total number of reads mapping
• Total RPK

First, take a quick look at the first five lines of the abundance.tsv for MT1.

There are five columns which give us information about the transcript abundances for our MT1 sample.

First, look for your gene of interest, PCHAS_1402500. Run this as a loop to grep the information for all five samples.

Now you have the length (eff_length) and counts (est_counts) for PCHAS_1402500 for each of your samples. Next, you need to get the total number of reads mapped to each of your samples. You can use a loop to do this.

In the loop below samtools flagstat gives you the number of mapped paired reads (reads with itself and mate mapped) and those where one read mapped but its mate didn't (singletons). If then uses grep to get the relevant lines and awk to add the mapped paired and singleton read totals together.

Finally, to calculate the TPM values, you need the total RPK for each of your samples. Again we use a loop. Notice the use of NR>2 in the awk command which tells it to skip the two header lines at the start of the file. You will also notice that we divide the eff_length by 1,000 so that it's in kilobases.

Q1: Using the abundance.tsv files generated by Kallisto and the information above, calculate the RPKM for PCHAS_1402500 in each of our five samples.

Q2: Using the abundance.tsv files generated by Kallisto and the information above, calculate the TPM for PCHAS_1402500 in each of our five samples.

Q3: Do these match the TPM values from Kallisto?

Yes.

Well, almost. They may be a couple out because we rounded up to make the calculations easier.

Q4: Do you think PCHAS_1402500 is differentially expressed between the MT and SBP samples?

Probably not.

If we were to look at only the counts and RPKM values then it appears there is an 8 fold difference between the MT and SBP samples. However, when we look at the TPM values, they are much closer and so differential expression is less likely.