# MA570 - Mapping and transcriptome assembly # The aim of this practical is to map RNA-Seq data to the human genome using TopHat and then to assemble the transcriptome using Cufflinks. # The data for this practical came from the ENCODE project which has RNA-Seq data for a number of cell types. # Today we will be using data from two cell types: # - H1hesc: Human embyronic stem cells. # - K562: Leukemia cancer cells. # The data we will be using is from human chromosome 18. # This is paired end data, with the left pair contained in files ending with _1.fq and the right in files ending in _2.fq # First, change directories to your own folder: cd #go to your own directory to begin the mapping ########################################################## # # # Step 1: Mapping reads to the human genome using TopHat # # # ########################################################## # These commands will map data onto the human genome (hg19). # The output will be placed in seperate folders for each run, into a file calles 'accepted_hits.bam' tophat -o H1hesc.chr18 --mate-inner-dist 50 --mate-std-dev 20 -p 5 -G /home/nextgen/data/RefSeq.hg19.chr18.23Feb.2013.gtf /home/nextgen/data/indexes/chr18 /home/nextgen/data/RNA_Seq/chr18/H1hesc.chr18_1.fq /home/nextgen/data/RNA_Seq/chr18/H1hesc.chr18_2.fq tophat -o K562.chr18 --mate-inner-dist 50 --mate-std-dev 20 -p 5 -G /home/nextgen/data/RefSeq.hg19.chr18.23Feb.2013.gtf /home/nextgen/data/indexes/chr18 /home/nextgen/data/RNA_Seq/chr18/K562.chr18_1.fq /home/nextgen/data/RNA_Seq/chr18/K562.chr18_2.fq # NOTE: TopHat may take a while to run on each dataset. # In the meantime, skip to Part 2, using pre-mapped data I have in /home/nextgen/peter # TopHat can be run with a nuimber of options. The options we used here are explained below: # --mate-inner-dist This is the expected (mean) inner distance between pairs. # For, example, for paired end runs with fragments selected at 200bp, where each end is 75bp, you should set this option to be 50 (as above). # --mate-std-dev The standard deviation for the distribution on inner distances between pairs. # -p This is a multi-threading options. Here we used 5 threads. # -G This specified a reference annotation file. # The file used here was downloaded from the UCSC genome browser and contains transcript annotations from the RefSeq database. ########################################### # # # Step 2: Quality control on mapped reads # # # ########################################### # Samtools is a widely used tool for BAM/SAM manipulations. It has many functions for processing/examining alignments. # Below are some examples: # Remove potential PCR duplicates samtools rmdup /home/nextgen/peter/H1hesc.chr18/accepted_hits.bam H1hesc.chr18.bam samtools rmdup /home/nextgen/peter/K562.chr18/accepted_hits.bam K562.chr18.bam # Check how many reads have aligned samtools view -c H1hesc.chr18.bam samtools view -c K562.chr18.bam # Some quality stats based on flags in the .BAM file samtools flagstat H1hesc.chr18.bam samtools flagstat K562.chr18.bam # More information about samtools can be found at http://samtools.sourceforge.net/samtools.shtml ################################################# # # # Step 3: Transcriptome assembly with Cufflinks # # # ################################################# cufflinks --multi-read-correct -p 10 -o H1hesc.chr18_cuffout ./H1hesc.chr18.bam cufflinks --multi-read-correct -p 10 -o K562.chr18_cuffout ./K562.chr18.bam # Options: # --milti-read-correct: Here cufflinks tries to deal with reads which map to multiple places in the genome. # Individual reads will sometimes be mapped to multiple positions in the genome due to sequence repeats and homology. # By default, Cufflinks will uniformly divide each multi-mapped read to all of the positions it maps to. # In other words, a read mapping to 10 positions will count as 10% of a read at each position. # With this option, cufflinks tries to 'rescue' such reads by re-estimating the abundances, dividing each multi-mapped read probabalistically # based on the initial abundance estimation of the genes it maps to and the inferred fragment length. ########################################## # # # Step 4: Merge the resulting assemblies # # # ########################################## ls -1 *cuffout/transcripts.gtf > assemblies.txt # This line creats a file assemblies.txt which has the path to the transcripts.gtf file we generated earlier. Cuffmerge will use this list to find these files. cuffmerge -s /home/nextgen/data/indexes/chr18.fa -p 10 -g /home/nextgen/data/RefSeq.hg19.chr18.23Feb.2013.gtf assemblies.txt # The -s option points to a copy of the genome, which is stored in fastA format. ################################################## # # # Step 5: Finding differentially expressed genes # # # ################################################## cuffdiff --multi-read-correct -p 10 -o cuffdiff_out merged_asm/merged.gtf ./H1hesc.chr18.bam ./K562.chr18.bam less cuffdiff_out/gene_exp.diff #-> to see the differentially expressed genes. ################################################## # # # Step 6: Have a look at the alignment using IGV # # # ################################################## # To have a look at the .bam file, we will first need to index the .bam file. # Samtools has a index function to do this. samtools index H1hesc.chr18.bam samtools index K562.chr18.bam # This will create .bam.bai files, which are the bam indexes. # Next, open IGV by typing igv # And then load the .bam file (not the .bam.bai file... IGV will automatically look for this)