Finally to get the counts

Now with all the bam files, indexed bam files and the bed file made previously we can get read counts per transcript. The moment you are anticipating is not far off. Make sure the bam.bai files are all in the same working directory.

We use BEDtools v2.24.0, with the multiBamCov option.

#PBS -P project
#PBS -N Syzygium_counts
#PBS -l nodes=1:ppn=20
#PBS -l walltime=01:00:00
#PBS -l pmem=24gb
#PBS -e ./Syzygium_counts.txt
#PBS -M email@sydney.edu.au
#PBS -m abe

# Load modules
module load bedtools

# Working directory
cd /working directory location

# Run bedtools script
multiBamCov -q 30 -p -bams AP0_coordSorted.bam AP1_coordSorted.bam AP2_coordSorted.bam -bed Syzygium.BED > Syzygium_counts.txt

The output from this will be a large text file with the numbers of reads that aligned to each transcript.

TRINITY_DN4630_c0_g1_i1 0 2449 68 0 18 30
TRINITY_DN4630_c0_g2_i1 0 457 2 0 4 4 0
TRINITY_DN4679_c0_g1_i1 0 410 2 0 2 0 0
TRINITY_DN4674_c0_g1_i1 0 227 4 0 0 0 0
TRINITY_DN4609_c0_g1_i1 0 265 0 0 0 0 0
TRINITY_DN4651_c0_g1_i1 0 1030 8 0 4 4

which you can put into a spreadsheet and name the columns. And the columns are in the order that the bam files are placed in the script. Of course you would include all the bam files for all samples and all times in the one run using bedtools.

Gene ID start end AP0 AP1 AP2
TRINITY_DN4630_c0_g1_i1 0 2449 68 0 18
TRINITY_DN4630_c0_g2_i1 0 457 2 0 4
TRINITY_DN4679_c0_g1_i1 0 410 2 0 2
TRINITY_DN4674_c0_g1_i1 0 227 4 0 0
TRINITY_DN4609_c0_g1_i1 0 265 0 0 0
TRINITY_DN4651_c0_g1_i1 0 1030 8 0 4

Aligning the reads to the transcriptome

Armed with the clustered fasta file and the trimmed raw reads you can now determine the read hits per transcript.

Alignment functions are supported within the Trinity platform, thereby incorporating the use of these third party tools, providing they are installed. Installed software required include; Bowtie (I used v1.1.2) and SAMtools (v1.2) to give outputs of bam files. Also required is Perl and of course Trinity – if running within the streamlined utilities.

#PBS -P project name
#PBS -N AP0_alignment
#PBS -l nodes=1:ppn=20
#PBS -l walltime=20:00:00
#PBS -l pmem=24gb
#PBS -e ./AP0_alignment.txt
#PBS -M email@sydney.edu.au
#PBS -m abe

# Load modules
module load perl
module load bowtie
module load java
module load samtools/1.2
module load trinity/2.1.1

# Working directory
cd /path to working directory

# Run bowtie script
/usr/local/trinity/2.1.1/util/bowtie_PE_separate_then_join.pl –seqType fq \
–left AP0_R1_pairedwithunpaired.trim.fq –right AP0_R2_pairedwithunpaired.trim.fq \
–target Syzygium.fasta –aligner bowtie \
— -p 4 –all –best –strata -m 300

From this will be output a folder called bowtie containing bam files and indexed bams. For the next step you will use the bam files from all samples and all times to determine the counts per transcript.I rename the bam files according the plant/time eg. AP0_coordSorted.bam and AP0_coordSorted.bam.bai

The incredible power of computing and what it means for biology

I want to express some of my wonder at the extraordinary power that is so readily available to us to examine highly complex biological interactions.

I find it incredible that we can manipulate and interpret the massive data-sets gene expression studies create. This can be done with ease using modern high performance computers. The point is that no-one could have done the sort of biology that I am doing just 10 or 15 years ago.

The large amounts of raw data from sequencing 8 plants at 3 moments in time (24 samples) came to around 48 Gigabytes of data. This is an enormous amount of biological information and each data-set is just a snap shot in time!!

Computers were co-opted into biology fairly recently and this has led to a rapid expansion in the possibilities available to understand life processes. One of the unexpected, but truly fantastic, side-effects of this merging of the sciences is the scrutiny of biology from a new computational logic perspective. It is not that biologists lack rigor but I think that the personality types attracted to biology have expanded – and this will lead to new approaches and more rapid advances in biological knowledge.

 

Running biopython scripts

As I am learning a bit about scripts and how to run them I have made recent errors which probably appear very simple to someone from a computing background.

The word I have learnt is ‘parse‘ which means:

  • to analyze (a sentence) into its parts and describe their syntactic roles.
  • analyze (a string or text) into logical syntactic components, typically in order to test conformability to a logical grammar.
  • examine or analyze minutely.

Another word is to ‘pass‘. You pass some instruction to the software to parse (analyze).

To run software you need:

  • the software and dependencies (other software) to be installed
  • script to tell the program what you want it to do
  • input files

I managed to find this very useful biopython script to make bed files from my newly made transcriptomes. Bed files are tab  delimited files that include the sequence ID, and the start and end location. For new transcriptomes the start is always 0 and the end is the length of the contig. These files can be used to gather the count data from your alignment files per sequence ID (gene).

http://ged.msu.edu/angus/tutorials-2013/rnaseq_bwa_counting.html

In case the website above is unavailable, here is the biopython script, slightly altered for the later version of Biopython which requires brackets around the print function:

import sys
from Bio import SeqIO

fasta_handle = open(sys.argv[1], "rU") #Open the fasta file specified by the user on the command line

#Go through all the records in the fasta file
for record in SeqIO.parse(fasta_handle, "fasta"):
    cur_id = record.id #Name of the current feature
    cur_length = len(record.seq) #Size of the current feature
    print (("%s\t%d\t%d") % (cur_id, 0, cur_length)) #Output the name, start, and end coordinates to the screen

This script will take any fasta file you give it and output the bed file. The important thing is to tell python.exe what to do with the input file. You need to pass the script to python.exe so it can parse and act on it. Therefore the order of your command must be correct – the .py script comes directly after the python.exe.

I placed the input file first and could not understand why I kept getting the error that my fasta file had syntax errors.

So in the command line window – I just run this in Windows cmd.exe by holding shift and right-clicking on the folder I want as my working directory;

Open the command window in the directory where your fasta file is, then input the directory path to Python34\python.exe then input the directory path to the script and output to your working directory.

So it will look something like this:

D:\Syzygium\>C:Python34\python.exe D:\scripts\make_bed_from_fasta.py Syzygium.fasta > Syzygium.bed

The output will look like this showing the length of each ‘gene’.

TRINITY_DN4638_c0_g1_i1 0 486
TRINITY_DN4638_c0_g2_i1 0 486
TRINITY_DN4608_c0_g1_i1 0 442
TRINITY_DN4687_c0_g1_i1 0 207
TRINITY_DN4630_c0_g1_i1 0 2449
TRINITY_DN4630_c0_g2_i1 0 457
TRINITY_DN4679_c0_g1_i1 0 410
TRINITY_DN4674_c0_g1_i1 0 227
TRINITY_DN4609_c0_g1_i1 0 265
TRINITY_DN4651_c0_g1_i1 0 1030
TRINITY_DN4626_c0_g1_i1 0 468
TRINITY_DN4626_c1_g1_i1 0 286
TRINITY_DN4653_c0_g2_i1 0 555
TRINITY_DN4653_c1_g1_i1 0 793
TRINITY_DN4648_c0_g1_i1 0 605
TRINITY_DN4632_c0_g1_i1 0 2026
TRINITY_DN4672_c0_g1_i1 0 239

 

Clustering the transcriptomes

I have now made eight transcriptomes. Four from the resistant plants and four from the susceptible. To be able to compare expression from different times, and between the different biological replicates, I have to make a single transcriptome for the plant species, in this case, Syzygium (Lilly Pilly). This is a bit concerning as my plants are grown from wild sourced seeds and are likely to be highly heterozygous.

My approach therefore was to use the CD-HIT-EST-2D (http://weizhongli-lab.org/cd-hit/) software to combine transcriptomes, initially at 0.95  similarity. I clustered the transcriptomes in pairs from the largest to the smallest until I had a single fasta file.

#PBS -P (project)
#PBS -N SYZ1_CD-HIT
#PBS -l nodes=1:ppn=20
#PBS -l walltime=00:15:00
#PBS -l pmem=24gb
#PBS -e ./SYZ1.error
#PBS -M (email address)
#PBS -m abe

# Load modules
module load cd-hit/4.6.1

# Working directory
Input path to directory

# Run CD-Hit script
cd-hit-est-2d -i BU.fasta -i2 BS.fasta -o SYZ1.fasta -c 0.95 -n 10 -d 0 -M 0 -T 0

I then sought out highly conserved genes using local BLAST in BioEdit. I used a couple of Eucalyptus grandis chitinase genes that I know to be well conserved from previous work, currently under review.

I found that the genes were present and intact in the individual transcriptomes but not in the combined one. I therefore ran the process with CD-HIT again with a higher stringency of 0.98.

This time the genes were present in the final transcriptome making me feel more confident about the alignment process.