##############################################################
#                                                            #
# MA570 - Genome Assembly and De Novo Transcriptome Assembly #
#                                                            #
##############################################################


###################
#                 #
# Genome assembly #
#                 #
###################

# In this practical we will look at a small genome assembly by attempting to assemble reads from human chromosome 22.
# For this we will be using single ended reads and the Abyss assembler.

# A parameter for this assmbly is kmer size. I've already determined the best kmer size for this to be 32bp

/home/sofia/software/bin/ABYSS -k 32 /home/nextgen/data/assembly/Abyss/NA12878.Chr22.fastq -o Chr22.assembled.fasta 

# Once assembled, have a look at the assembly to inspect the contigs.

less Chr22.assembly.fasta

# As you can see, we get many short contigs. Try calculating the N50 for this assembly.

perl /cm/shared/apps/Trinity/r2012-10-05/util/N50.pl Chr22.assembled.fasta


##################################
#                                #
# De Novo Transcriptome Assembly #
#                                #
##################################

# In this practical we will be using paired end RNA-Seq reads to assemble a transcriptome from scratch.
# That is, without a reference genome.

# To do this we will be running the Trinity assembler.
# Trinity consists of three main steps; Inchworm, Chrysalis and Butterfly.
# What each one does is below. 

# Inchworm assembles the RNA-seq data into the unique sequences of transcripts, 
# often generating full-length transcripts for a dominant isoform, 
# but then reports just the unique portions of alternatively spliced transcripts.

# Chrysalis clusters the Inchworm contigs into clusters and constructs complete de Bruijn graphs for each cluster. 
# Each cluster represents the full transcriptonal complexity for a given gene (or sets of genes that share sequences in common). 
# Chrysalis then partitions the full read set among these disjoint graphs. 

# Butterfly then processes the individual graphs in parallel, 
# tracing the paths that reads and pairs of reads take within the graph, 
# ultimately reporting full-length transcripts for alternatively spliced isoforms, 
# and teasing apart transcripts that corresponds to paralogous genes. 

# Run Trinity as follows:

Trinity.pl --seqType fq --JM 10G --left /home/nextgen/data/assembly/Trinity/Chr18.left.fastq --right /home/nextgen/data/assembly/Trinity/Chr18.left.fastq --CPU 6

# Options:
# --seqtype fq : Tells Trinity the data is in fastQ format
# --JM 10G : The total memory given to Trinity for kmer searching (10 gigabytes here)
# --left : left read pairs
# --right : right read pairs
# --CPU : Number of processors to use


#################################
#                               #
# Examine the output of Trinity #
#                               #
#################################

# The output for Trinity is in the directory trinity_out_dir/
# Go into this folder 

cd trinity_out_dir/

# The assembled transcripts are in a file Trinity.fasta
# Count the number of transcripts assembled.

grep -c ">" Trinity.fasta

# Have a look at the transcripts.

less Trinity.fasta

# Calculate the N50 for this assembly.

perl /cm/shared/apps/Trinity/r2012-10-05/util/N50.pl Trinity.fasta


####################################################
#                                                  #
# Map reads back to the transcriptome and quantify #
#                                                  #
####################################################

# The next step involves mapping reads back to the newly assembled transcriptome using Bowtie
# Trinity has a script that will do that for us.

perl /cm/shared/apps/Trinity/r2012-10-05/util/alignReads.pl --left /home/nextgen/data/assembly/Trinity/Chr18.left.fastq  --right /home/nextgen/data/assembly/Trinity/Chr18.right.fastq --seqType fq --target Trinity.fasta --aligner bowtie

# Once this has beed completed, run another of Trinitys scripts to estimate abundance of each transcript.

cd bowtie_out/
perl /cm/shared/apps/Trinity/r2012-10-05/util/RSEM_util/run_RSEM.pl --transcripts /home/nextgen/peter/Trinity/trinity_out_dir/Trinity.fasta --name_sorted_bam bowtie_out.nameSorted.bam --paired

# This generates a number of files, including RSEM.genes.results, RSEM.isoforms.results which contain expression values in FPKM.
# This can be used for downstream analysis.