\n
Summary<\/h1>\n
<\/p>\n
- \n
- Preambule : 0. Going to the i-Trop cluster -
ssh,srun,scp<\/code><\/a><\/li>\n- Practice 1:# 1. Check Reads Quality<\/a><\/li>\n
- Practice 2:# 2. Performing a de novo RNA-Seq assembly with
trinity<\/code><\/a><\/li>\n- Practice 3:# 3. Assessing transcriptome assembly quality<\/a><\/li>\n
- Practice 4:# 4. Differential Expression Analysis (DE)<\/a><\/li>\n
- License<\/a><\/li>\n<\/ul>\n
\n<\/a><\/p>\n
0. Going to the i-Trop cluster -
ssh,srun,scp<\/code><\/h1>\nDataset used in this practical comes from<\/p>\n
- \n
- ref : https:\/\/www.ncbi.nlm.nih.gov\/pmc\/articles\/PMC3488244\/<\/a><\/li>\n
- data : NCBI SRA database under accession number SRS307298 S. cerevisiae<\/li>\n
- Genome size of S. cerivisiae : 12M (12.157.105) (https:\/\/www.yeastgenome.org\/strain\/S288C#genome_sequence<\/a>)<\/li>\n<\/ul>\n
In this session, we will analyze RNA-seq data from one sample of S. cerevisiae (NCBI SRA
\nSRS307298). It is from two different origin (CENPK and Batch), with three biological replications for each
\norigin (rep1, rep2 and rep3).<\/p>\nConnection to the i-Trop Cluster through
ssh<\/code> mode<\/h3>\nWe will work on the i-Trop Cluster with a "supermem" node using SLURM scheduler.<\/p>\n
ssh formationX@bioinfo-master.ird.fr<\/code><\/pre>\nOpening an interactive bash session on the node25 (supermem partition) -
srun -p partition --pty bash -i<\/code><\/h3>\nRead this survival document containig basic commands to SLURM (https:\/\/southgreenplatform.github.io\/trainings\/slurm\/<\/a>)<\/p>\n
srun -p supermem --mem 50G --time 20:00:00 --cpus-per-task 2 --pty bash -i<\/code><\/pre>\nPrepare input files<\/h3>\n
- \n
- Create your subdirectory in the scratch file system \/scratch. In the following, please replace X with your own user ID number in formationX.<\/li>\n<\/ul>\n
cd \/scratch\nmkdir formationX\ncd formationX<\/code><\/pre>\n- \n
- Copy the exercise files from the shared location to your scratch directory (it is essential that all
\ncalculations take place here)<\/li>\n<\/ul>\nscp -r nas:\/data2\/formation\/TP-trinity\/SRA_SRS307298\/RAWDATA\/ \/scratch\/formationX\/<\/code><\/pre>\n- \n
- When the files transfer is finished, verify by listing the content of the current directory and the subdirectory RAWDATA with the
\ncommandls -al<\/code>. You should see 12 gzipped read files in a listing, thesamples.txt<\/code> file and therun_trinity.sh<\/code> bash script. <\/li>\n<\/ul>\n[orjuela@node25 RAWDATA]$ more samples.txt \nCENPK CENPK_rep1 PATH\/SRR453569_1.fastq.gz PATH\/SRR453569_2.fastq.gz\nCENPK CENPK_rep2 PATH\/SRR453570_1.fastq.gz PATH\/SRR453570_2.fastq.gz\nCENPK CENPK_rep3 PATH\/SRR453571_1.fastq.gz PATH\/SRR453571_2.fastq.gz\nBatch Batch_rep1 PATH\/SRR453566_1.fastq.gz PATH\/SRR453566_2.fastq.gz\nBatch Batch_rep2 PATH\/SRR453567_1.fastq.gz PATH\/SRR453567_2.fastq.gz\nBatch Batch_rep3 PATH\/SRR453568_1.fastq.gz PATH\/SRR453568_2.fastq.gz<\/code><\/pre>\n
\n<\/a><\/p>\n
1. Check Reads Quality<\/h1>\n
FastQC perform some simple quality control checks to ensure that the raw data looks good and there are no problems or biases in data which may affect how user can usefully use it. http:\/\/www.bioinformatics.babraham.ac.uk\/projects\/fastqc\/<\/a><\/p>\n
# make a fastqc repertoty\ncd ..\/\nmkdir FASTQC; cd FASTQC\n\n#charge modules \nmodule load bioinfo\/FastQC\/0.10.1\n\n# run fastqc in the whole of samples\nfastqc -t 2 \/scratch\/formationX\/RAWDATA\/*.gz -o \/scratch\/formationX\/FASTQC\/<\/code><\/pre>\nMultiqc is a modular tool to aggregate results from bioinformatics analyses across many samples into a single report. Use this tool to visualise results of quality. https:\/\/multiqc.info\/<\/a><\/p>\n
#charge modules \nmodule load bioinfo\/multiqc\/1.7\n\n#launch Multiqc to create a report in html containing the whole of informations generated by FastQC\nmultiqc .\n\n#transfer results to your cluster home \nscp -r multiqc* nas:\/home\/formationX\/\n\n# transfert results to your local machine by scp or filezilla\nscp -r formationX@bioinfo-nas.ird.fr:\/home\/formationX\/multiqc* .\/\n\n# open in your favorite web navigator\nfirefox multiqc_report.html .<\/code><\/pre>\nIn this practice, reads quality is ok. You need to observe sequences and check biases. To remove adaptors and primers you can use Trimmomatic. Use PRINSEQ2 to detect Poly A\/T tails and low complexity reads. Remove contaminations with SortMeRNA, riboPicker or DeconSeq.<\/p>\n
\n<\/a><\/p>\n
2. Performing a de novo RNA-Seq assembly with
trinity<\/code><\/h1>\n2.1 Running Trinity with trimmomatic and reads normalisation<\/h2>\n
Preparing assembly sample file and check parameters of trinity assembler<\/h4>\n
Observe
run_trinity.sh<\/code> script and adapt to your formation number. This script is not sending in sbatch mode because it could be take time in this training. If you need launch it in a slurm cluster you can use this version containing SBATCH commands run_trinity.slurm<\/a> or adapt it to SGE.<\/p>\n[orjuela@node25 RAWDATA]$ more run_trinity.sh \n\n# loading modules\nmodule load bioinfo\/samtools\/1.9\nmodule load bioinfo\/trinityrnaseq\/2.8.5\nmodule load bioinfo\/Trimmomatic\/0.33\n\n# changing PATH to current directory in samples file\nsed -i 's|PATH|'$PWD'|ig' samples.txt \n\n# Running trinity assembly\n# Trinity --seqType fq --max_memory 50G --CPU 2 --samples_file samples.txt --output ..\/TRINITY_OUT \n\n# Running Trinity with trimmomatic and normalisation\nTrinity --seqType fq --max_memory 50G --CPU 2 --trimmomatic --quality_trimming_params 'ILLUMINACLIP:\/usr\/local\/Trimmomatic-0.33\/adapters\/TruSeq2-PE.fa:2:30:10 ILLUMINACLIP:\/scratch\/formationX\/RAWDATA\/adapt-125pbLib.txt:2:30:10 SLIDINGWINDOW:5:20 LEADING:5 TRAILING:5 MINLEN:25 HEADCROP:10' --normalize_by_read_set --samples_file samples.txt --output ..\/TRINITY_OUT<\/code><\/pre>\nRunning assembly :<\/p>\n
bash run_trinity.sh > ..\/trinity.log & <\/code><\/pre>\nAll screen output (info messages and error messages, if any) will be saved in the file
\n..\/trinity.log<\/code>. The script will start executing in the background (the & at the end), so
\nthat the terminal will return to the prompt right after you hit \u201cEnter\u201d. <\/p>\nYou can use
jobs<\/code> to monitor jobs but if you logout the program does not keep running.<\/p>\nWhile running you can examine steps
tail -f ..\/trinity.log<\/code>.<\/p>\nAfter trimmomatic and reads normalisation, Three stages are done by Trinity<\/p>\n
- \n
- Jellyfish: Extracts and counts K-mers (K=25) from reads<\/li>\n
- Inchworm: Assembles initial contigs by \u201cgreedily\u201d extending sequences with most abundant K-mers<\/li>\n
- Chrysalis: Clusters overlapping Inchworm contigs, builds de Bruijn graphs for each cluster, partitions reads between clusters<\/li>\n
- Butterfly: resolves alternatively spliced and paralogous transcripts independently for each cluster (in parallel)<\/li>\n<\/ol>\n
Ending up assembly and downloading assembly results<\/h4>\n
- \n
- \n
WARNING !: This job is expected to run 12 hours. Kill your job using
fg<\/code> andctl+c<\/code> <\/p>\n<\/li>\n<\/ul>\njobs = shows jobs running\n fg = foreground\n bf = background\n ctrl+z = send job to background\n ctrl+c = kill job in ground<\/code><\/pre>\n- \n
- Recover assembly results (only Trinity.fasta) generated by trinity from the shared location to your scratch directory : <\/li>\n<\/ul>\n
scp -r nas:\/data2\/formation\/TP-trinity\/TRINITY_OUT\/TRINITY_OUT\/Trinity.fasta \/scratch\/formationX\/\n\ncd \/scratch\/formationX\/TRINITY_OUT<\/code><\/pre>\n- \n
- Upon successful completion of Trinity, the assembled transcriptome is written to the FASTA file called
\nTrinity.fasta<\/code> located in the output directory\/scratch\/formationX\/TRINITY_OUT<\/code>.<\/li>\n<\/ul>\n
\n<\/a><\/p>\n
3. Assessing transcriptome assembly quality<\/h2>\n
3.1 Getting basic Assembly metrics with the trinity script
TrinityStats.pl<\/code><\/h2>\nRunning trinity script<\/h5>\n
Trinity.fasta contains transcripts to be evaluated, annotated, and used in downstream analysis of
\nexpression. In this exercise, we only concentrate on basis statistics of the assembled transcriptome,
\nwhich can be obtained using a Trinity utility script TrinityStats.pl. <\/p>\n#declare the path_to_trinity variable\npath_to_trinity=\/usr\/local\/trinityrnaseq-2.8.5\/\n#check path\necho $path_to_trinity\n#launch assembly metrics script\n$path_to_trinity\/util\/TrinityStats.pl \/scratch\/formationX\/TRINITY_OUT\/Trinity.fasta > trinityStats.log<\/code><\/pre>\nParsing the results generated by the script<\/h5>\n
The output file generated (trinityStats.log) will contain basic information about contig length distributions, based on all transcripts and only on the longest isoform per gene. <\/p>\n
Besides average and median contig lengths, also given are quantities N10 through N50. Nx is the smallest contig length such that (x\/100)% of all assembled bases are in contigs longer than Nx. Specifically, N50 is the contig length such that half of all assembly sequence is contained in contigs longer than that. In whole genome assembly, N50 is often used as a measure (one of many) of assembly quality, since the longer the contigs, the better the assembly. In the case of transcriptome, contig lengths should be correct, which does not imply \u201clarge\u201d. If it falls in the right ballpark (about 1000-1,500), N50 can still be used as a check on overall \u201csanity\u201d of the transcriptome assembly.<\/p>\n
[orjuela@node25 orjuela]$ more trinityStats.log \n\n################################\n## Counts of transcripts, etc.\n################################\nTotal trinity 'genes': 7600\nTotal trinity transcripts: 8616\nPercent GC: 39.37\n\n########################################\nStats based on ALL transcript contigs:\n########################################\n\n Contig N10: 11180\n Contig N20: 8394\n Contig N30: 6793\n Contig N40: 5553\n Contig N50: 4736\n\n Median contig length: 506\n Average contig: 1909.88\n Total assembled bases: 16455500\n\n#####################################################\n## Stats based on ONLY LONGEST ISOFORM per 'GENE':\n#####################################################\n\n Contig N10: 9989\n Contig N20: 7558\n Contig N30: 6192\n Contig N40: 5131\n Contig N50: 4194\n\n Median contig length: 423\n Average contig: 1625.41\n Total assembled bases: 12353135\n<\/code><\/pre>\n3.2 Reads mapping back rate and abundance estimation using the trinity script
align_and_estimate_abundance.pl<\/code><\/h2>\nRead congruency is an important measure in determining assembly accuracy. Clusters of read pairs that align incorrectly are strong indicators of mis-assembly. A typical "good" assembly has ~80% reads mapping to the assembly and ~80% are properly paired.<\/p>\n
Several tools can be used to calculate reads mapping back rate over
Trinity.fasta<\/code> assembly : bwa, bowtie2 (mapping), kallisto, salmon (pseudo mapping). Quantify read counts for each gene\/isoform can be calculate. Mapping and quantification can be obtained by using the --est_method argument into thealign_and_estimate_abundance.pl<\/code> script.<\/p>\nWe will performing this analyses step successively with the
align_and_estimate_abundance.pl<\/code> script :<\/p>\n- \n
- Pseudomapping methods (kallisto or salmon) are faster than mapping based. So firstly we will use salmon to pseudoalign reads from sample to the reference and quantify abondance. <\/li>\n
- Then, we will use bowtie2 and RSEM to align and quantify read counts for each gene\/isoform.<\/li>\n<\/ul>\n
Preparing data and environment for our following analysis<\/h5>\n
#create a directory dans \/scratch\/formationX\/\nmkdir ALIGN_AND_ABUNDANCE\ncd ALIGN_AND_ABUNDANCE\n\n# Loading modules\nmodule load system\/perl\/5.24.0\nmodule load bioinfo\/trinityrnaseq\/2.8.5\nmodule load bioinfo\/bowtie2\/2.2.9\nmodule load bioinfo\/express\/1.5.1\nmodule load bioinfo\/kallisto\/0.43.1\nmodule load bioinfo\/RSEM\/1.0\nmodule load bioinfo\/salmon\/0.10.2\nmodule load bioinfo\/samtools\/1.7\n\n# modifying samples file path because quality-trimmed your reads using the --trimmomatic parameter in Trinity\ncp \/scratch\/formationX\/RAWDATA\/samples.txt .\nsed -i 's|RAWDATA|TRINITY_OUT|ig' samples.txt\nsed -i 's|.fastq.gz|.fastq.gz.P.qtrim.gz|ig' samples.txt\n\n# define variables\npath_to_trinity=\/usr\/local\/trinityrnaseq-2.8.5\/\nfasta=\/scratch\/formationX\/TRINITY_OUT\/Trinity.fasta\nsamplesfile=\/scratch\/formationX\/ALIGN_AND_ABUNDANCE\/samples.txt #modified samples file<\/code><\/pre>\nLaunch Salmon<\/h5>\n
# create a salmon_outdir and go on\nmkdir salmon_outdir; cd salmon_outdir\n\n# salmon\nperl $path_to_trinity\/util\/align_and_estimate_abundance.pl \\\n--transcripts $fasta \\\n--seqType fq \\\n--samples_file $samplesfile \\\n--est_method salmon \\\n--trinity_mode \\\n--prep_reference > salmon_align_and_estimate_abundance.log 2>&1 &<\/code><\/pre>\nCheck reads percentage mapped to Trinity.fasta reference by sample for salmon and bowtie-rsem est method.<\/h4>\n
#salmon results\n[orjuela@node25 salmon_outdir]$ grep 'Mapping rate =' *\/logs\/salmon_quant.log \nBatch_rep1\/logs\/salmon_quant.log:[2019-09-17 15:54:32.923] [jointLog] [info] Mapping rate = 96.3027%\nBatch_rep2\/logs\/salmon_quant.log:[2019-09-17 15:55:09.922] [jointLog] [info] Mapping rate = 96.3992%\nBatch_rep3\/logs\/salmon_quant.log:[2019-09-17 15:55:39.246] [jointLog] [info] Mapping rate = 96.7795%\nCENPK_rep1\/logs\/salmon_quant.log:[2019-09-17 15:52:56.609] [jointLog] [info] Mapping rate = 94.4516%\nCENPK_rep2\/logs\/salmon_quant.log:[2019-09-17 15:53:30.146] [jointLog] [info] Mapping rate = 95.0558%\nCENPK_rep3\/logs\/salmon_quant.log:[2019-09-17 15:54:00.196] [jointLog] [info] Mapping rate = 95.0639%<\/code><\/pre>\nOPTIONAL: Launch bowtie2-rsem<\/h5>\n
- \n
- \n
WARNING !: this job can take a lot ~1h30 by sample. This step will not run in this practice.<\/p>\n<\/li>\n<\/ul>\n
# create a bowtie2-rsem_outdir and go on\ncd \/scratch\/formationX\/ALIGN_AND_ABUNDANCE\/\nmkdir bowtie2-rsem_outdir; cd bowtie2-rsem_outdir\n\n# runnign align_and_estimate_abundance in bowtie2-rsem mode\nperl $path_to_trinity\/util\/align_and_estimate_abundance.pl \\\n--transcripts $fasta \\\n--seqType fq \\\n--samples_file $samplesfile \\\n--est_method RSEM --aln_method bowtie2 \\\n--trinity_mode \\\n--prep_reference \\\n--coordsort_bam > bowtie-rsem_align_and_estimate_abundance.log 2>&1 &<\/code><\/pre>\nHere, we show you results with bowtie2-rsem OPTIONAL part<\/p>\n
#bowtie-rsem results\n[orjuela@node25 bowtie2-rsem_outdir]$ grep -B 6 "overall alignment rate" bowtie-rsem_align_and_estimate_abundance.log\nCMD: set -o pipefail && bowtie2 --no-mixed --no-discordant --gbar 1000 --end-to-end -k 200 -q -X 800 -x \/scratch\/orjuela\/TRINITY_OUT\/Trinity.fasta.bowtie2 -1 \/scratch\/orjuela\/TRINITY_OUT\/SRR453569_1.fastq.gz.P.qtrim.gz -2 \/scratch\/orjuela\/TRINITY_OUT\/SRR453569_2.fastq.gz.P.qtrim.gz -p 4 | samtools view -@ 4 -F 4 -S -b | samtools sort -@ 4 -n -o bowtie2.bam \n3666948 reads; of these:\n 3666948 (100.00%) were paired; of these:\n 497221 (13.56%) aligned concordantly 0 times\n 2058938 (56.15%) aligned concordantly exactly 1 time\n 1110789 (30.29%) aligned concordantly >1 times\n86.44% overall alignment rate\n--\n ...<\/code><\/pre>\n3.3. Expression matrix construction<\/h2>\n
Combine read count from all samples into a matrix, and normalize the read count using the TMM method. This command will take in RSEM output files from each sample, and combine them into a single matrix file.<\/p>\n
# go to salmon results repertory\ncd \/scratch\/formationX\/ALIGN_AND_ABUNDANCE\/salmon_outdir\n\n#load modules\nmodule load bioinfo\/R\/3.3.3\nmodule load system\/perl\/5.24.0\n\n#declare bash variables\npath_to_trinity=\/usr\/local\/trinityrnaseq-2.8.5\/\n\n# calculate expression matrix (if salmon use quant.sf files, if kallisto use abundance.tsv files)\n$path_to_trinity\/util\/abundance_estimates_to_matrix.pl \\\n--est_method salmon \\\n--out_prefix Trinity_trans \\\n--name_sample_by_basedir \\\n--gene_trans_map none \\\nCENPK_rep1\/quant.sf \\\nCENPK_rep2\/quant.sf \\\nCENPK_rep3\/quant.sf \\\nBatch_rep1\/quant.sf \\\nBatch_rep2\/quant.sf \\\nBatch_rep3\/quant.sf <\/code><\/pre>\nYou have to obtain two matrices:
\nThe firts one containing the estimated counts\u0300Trinity_trans.isoform.counts.matrix<\/code> and the second one containing the TPM expression values that are cross-sample normalized using the TMM methodTrinity_trans.TMM.EXPR.matrix<\/code>.
\nTMM normalization assumes that most transcripts are not differentially expressed, and linearly scales the expression values of samples to better enforce this property.<\/p>\nIn both files, each column represents a sample, and each row represents a gene, the values are either the raw read counts or normalized FPKM values. The \u201ccounts\u201d file will be used for differentially expressed gene identification, and the \u201cfpkm\u201d file will be used for clustering analysis. By default, the fpkm file is normalized with TMM method. <\/p>\n
[orjuela@node25 salmon_outdir]$ ll Trinity_trans.*\ntotal 1,4M\n-rw-r--r-- 1 orjuela borea-equipe7 539K 29 ao\u00fbt 16:41 Trinity_trans.isoform.TMM.EXPR.matrix\n-rw-r--r-- 1 orjuela borea-equipe7 386K 29 ao\u00fbt 16:41 Trinity_trans.isoform.counts.matrix\n-rw-r--r-- 1 orjuela borea-equipe7 474K 29 ao\u00fbt 16:41 Trinity_trans.isoform.TPM.not_cross_norm\n-rw-r--r-- 1 orjuela borea-equipe7 355 29 ao\u00fbt 16:41 Trinity_trans.isoform.TPM.not_cross_norm.TMM_info.txt\n-rw-r--r-- 1 orjuela borea-equipe7 542 29 ao\u00fbt 16:41 Trinity_trans.isoform.TPM.not_cross_norm.runTMM.R\n\n[orjuela@node25 salmon_outdir]$ head *.matrix\n==> Trinity_trans.isoform.counts.matrix <==\n CENPK_rep1 CENPK_rep2 CENPK_rep3 Batch_rep1 Batch_rep2 Batch_rep3\nTRINITY_DN3130_c0_g1_i1 0.000000 4.000000 0.000000 0.000000 0.000000 0.000000\nTRINITY_DN986_c0_g1_i4 190.000000 210.677452 319.000000 353.571187 420.080043 292.000000\nTRINITY_DN4062_c0_g2_i1 0.000000 13.000000 0.000000 0.000000 0.000000 0.000000\nTRINITY_DN4840_c0_g1_i1 0.000000 10.000000 0.000000 0.000000 0.000000 0.000000\nTRINITY_DN1624_c0_g1_i1 0.000000 5.000000 0.000000 0.000000 0.000000 0.000000\nTRINITY_DN3492_c0_g1_i1 0.000000 5.000000 0.000000 0.000000 0.000000 0.000000\nTRINITY_DN1386_c0_g1_i1 671.000000 724.000000 997.000000 1325.000000 1726.000000 1229.000000\nTRINITY_DN6688_c0_g1_i1 0.000000 5.000000 0.000000 0.000000 0.000000 0.000000\nTRINITY_DN1126_c0_g2_i1 29.000000 42.000000 51.000000 14.000000 7.000000 5.000000\n\n==> Trinity_trans.isoform.TMM.EXPR.matrix <==\n CENPK_rep1 CENPK_rep2 CENPK_rep3 Batch_rep1 Batch_rep2 Batch_rep3\nTRINITY_DN3130_c0_g1_i1 0.000 32.382 0.000 0.000 0.000 0.000\nTRINITY_DN986_c0_g1_i4 141.255 131.008 158.103 195.483 180.432 173.943\nTRINITY_DN4062_c0_g2_i1 0.000 20.291 0.000 0.000 0.000 0.000\nTRINITY_DN4840_c0_g1_i1 0.000 23.121 0.000 0.000 0.000 0.000\nTRINITY_DN1624_c0_g1_i1 0.000 22.187 0.000 0.000 0.000 0.000\nTRINITY_DN3492_c0_g1_i1 0.000 30.511 0.000 0.000 0.000 0.000\nTRINITY_DN1386_c0_g1_i1 96.472 86.724 94.337 140.804 141.542 140.796\nTRINITY_DN6688_c0_g1_i1 0.000 23.740 0.000 0.000 0.000 0.000\nTRINITY_DN1126_c0_g2_i1 9.819 11.863 11.412 3.511 1.357 1.351<\/code><\/pre>\n[orjuela@node25 salmon_outdir]$ more Trinity_trans.isoform.TPM.not_cross_norm.TMM_info.txt\ngroup lib.size norm.factors eff.lib.size\nCENPK_rep1 999995 1.14115425493455 1141148.54916327\nCENPK_rep2 999993 0.905154412329735 905148.076248848\nCENPK_rep3 1000003 1.0939694255262 1093972.70743448\nBatch_rep1 999992 0.989872299588587 989864.38061019\nBatch_rep2 999990 0.950965458077021 950955.94842244\nBatch_rep3 1000015 0.940121262909323 940135.364728266<\/code><\/pre>\n3.4 Compute N50 based on the top-most highly expressed transcripts (Ex50)<\/h2>\n
$path_to_trinity\/util\/misc\/contig_ExN50_statistic.pl Trinity_trans.isoform.TMM.EXPR.matrix \/scratch\/formationX\/TRINITY_OUT\/Trinity.fasta > ExN50.stats\n\n#Plotting ExN50\n$path_to_trinity\/util\/misc\/plot_ExN50_statistic.Rscript ExN50.stats\n\n#If you want to know, how many transcripts correspond to the Ex 90 peak, you could:\ncat Trinity_trans.isoform.TMM.EXPR.matrix.E-inputs | egrep -v ^\\# | awk '$1 <= 90' | wc -l\n1761\n\n# or 50\ncat Trinity_trans.isoform.TMM.EXPR.matrix.E-inputs | egrep -v ^\\# | awk '$1 <= 50' | wc -l\n163<\/code><\/pre>\nTRANSFERT: Observe plots. Remember, you have to transfert *.pdf files to your home before of transfering into your local machine.<\/p>\n
#transfering plots\ncp *.pdf \/home\/formationX\/\n# from your local machine\nscp formationX@bioinfo-nas.ird.fr:*.pdf .<\/code><\/pre>\n![\"\"](\"https:\/\/southgreenplatform.github.io\/trainings\/\/images\/Ex50.png\")
<\/p>\n3.5 Quantifying completness using
BUSCO<\/code><\/h2>\nAssessing gene space is a core aspect of knowing whether or not you have a good assembly.<\/p>\n
Benchmarking Universal Single-Copy Orthologs (BUSCO) sets are collections of orthologous groups with near-universally-distributed single-copy genes in each species, selected from OrthoDB root-level orthology delineations across arthropods, vertebrates, metazoans, fungi, and eukaryotes. BUSCO groups were selected from each major radiation of the species phylogeny requiring genes to be present as single-copy orthologs in at least 90% of the species; in others they may be lost or duplicated, and to ensure broad phyletic distribution they cannot all be missing from one sub-clade. The species that define each major radiation were selected to include the majority of OrthoDB species, excluding only those with unusually high numbers of missing or duplicated orthologs, while retaining representation from all major sub-clades. Their widespread presence means that any BUSCO can therefore be expected to be found as a single-copy ortholog in any newly-sequenced genome from the appropriate phylogenetic clade.<\/p>\n
Running Busco on Trinity.fasta assembly<\/h5>\n
#make a BUSCO directory\ncd \/scratch\/formationX\/\nmkdir BUSCO; cd BUSCO\n# download the eukaryota database from https:\/\/busco.ezlab.org\/\n# There is a lot of other databases usable for this dataset (ex: Fungi gene set or Saccaromyceta gene set, with more genes to retrieve, so longer to run)\n\nwget https:\/\/busco.ezlab.org\/datasets\/eukaryota_odb9.tar.gz\n# uncompress repertory\ntar zxvf eukaryota_odb9.tar.gz\n\n# define variables\nbusco=\/usr\/local\/BUSCO-3.0.2\/scripts\/run_BUSCO.py\npath_to_trinity=\/usr\/local\/trinityrnaseq-2.8.5\/\nfasta=\/scratch\/formationX\/TRINITY_OUT\/Trinity.fasta\nLINEAGE=\/scratch\/formationX\/BUSCO\/eukaryota_odb9\n\n# run busco\n\/usr\/local\/python-3.6.5\/bin\/python $busco -i $fasta -l $LINEAGE -m transcriptome -c 2 -o trinity_busco > busco.log &\n\n# check busco running using \njobs\n# or \ntail -f busco.log<\/code><\/pre>\nDisplaying the short summary table generated by BUSCO<\/h5>\n
[orjuela@node25 run_trinity_busco]$ more run_trinity_busco\/short_summary_trinity_busco.txt \n# BUSCO version is: 3.0.2 \n# The lineage dataset is: eukaryota_odb9 (Creation date: 2016-11-02, number of species: 100, number of BUSCOs: 303)\n# To reproduce this run: python \/usr\/local\/BUSCO-3.0.2\/scripts\/run_BUSCO.py -i \/scratch\/formation1\/TRINITY_OUT\/Trinity.fasta -o trinity_busco_euk -l \/scratch\/formation1\/BUSCO\/eukaryota_odb9\/ -m transcriptome -c 2\n\n#\n# Summarized benchmarking in BUSCO notation for file \/scratch\/orjuela\/TRINITY_OUT\/Trinity.fasta\n# BUSCO was run in mode: transcriptome\n\n C:93.4%[S:74.6%,D:18.8%],F:0.0%,M:6.6%,n:303\n\n 283 Complete BUSCOs (C)\n 226 Complete and single-copy BUSCOs (S)\n 57 Complete and duplicated BUSCOs (D)\n 0 Fragmented BUSCOs (F)\n 20 Missing BUSCOs (M)\n 303 Total BUSCO groups searched\n<\/code><\/pre>\n3.6 BLASTX comparison to known protein sequences database<\/h2>\n
Performing a blastx againt the swissprot database (the manually annotated and reviewed section of the UniProt Knowledgebase)<\/h5>\n
- \n
- \n
WARNING !: This step will not run in this practice because it will be inclued in annotation script !! But we give you commands lines to lauch with your data.<\/p>\n<\/li>\n<\/ul>\n
#load modules\nmodule load bioinfo\/blast\/2.4.0+\n\n#define variables\npath_to_blast=\/usr\/local\/ncbi-blast-2.4.0+\/bin\/\npath_to_trinity=\/usr\/local\/trinityrnaseq-2.8.5\/\n\n# create BLASTX repertory\nmkdir \/scratch\/formationX\/BLASTX\ncd \/scratch\/formationX\/BLASTX\n\n# First, we downloaded and indexed the database from ftp:\/\/ftp.ebi.ac.uk\/pub\/databases\/uniprot\/current_release\/knowledgebase\/complete\/uniprot_sprot.fasta.gz, \n# uniprotdatabase is available on cluster at (with the correct index files)\nuniprot=\/data\/projects\/banks\/Uniprot\/uniprot_sprot.fasta\n\n#Then, we ran BLASTX to get the top match hit:\n$path_to_blast\/blastx \\\n-db $uniprot \\\n-query $fasta \\\n-num_threads 2 \\\n-max_target_seqs 1 \\\n-outfmt 6 \\\n-evalue 1e-20 > SwissProt_1E20_Trinity_blastx.outfmt6 &\n\n#if you want to test blastX practice, please recover results from `\/data2\/formation\/TP-trinity\/TRINITY_OUT\/BLASTX\/`\n\n#Finally, we examined the percent of alignment coverage:\n$path_to_trinity\/util\/misc\/blast_outfmt6_group_segments.pl SwissProt_1E20_Trinity_blastx.outfmt6 $fasta $uniprot > SwissProt_1E20_Trinity_blastx.outfmt6.grouped\n\n$path_to_trinity\/util\/misc\/blast_outfmt6_group_segments.tophit_coverage.pl SwissProt_1E20_Trinity_blastx.outfmt6.grouped \\\n> SwissProt_1E20_Trinity_blastx.outfmt6.grouped.output<\/code><\/pre>\nAnalyzing blast results<\/h5>\n
[orjuela@node25 BLASTX]$ more SwissProt_1E20_Trinity_blastx.outfmt6.grouped.output\n#hit_pct_cov_bin count_in_bin >bin_below\n100 242 242\n90 8 250\n80 4 254\n70 5 259\n60 4 263\n50 5 268\n40 4 272\n30 2 274\n20 2 276\n10 1 277<\/code><\/pre>\n8616 transcripts from
Trinity.fasta<\/code> file were blasted against the swissprot database.<\/p>\n440 hits were reported in
SwissProt_1E20_Trinity_blastx.outfmt6<\/code> file.<\/p>\nIn
SwissProt_1E20_Trinity_blastx.outfmt6.grouped.output<\/code> file we can observed for example that 242 sequences were found with a 100% identity to an uniprot protein (count_in_bin). bin_below column represent a accumulative number of sequences.<\/p>\n
\n<\/a><\/p>\n
4. Differential Expression Analysis (DE)<\/h1>\n
4.1 Examine your data and your experimental replicates before DE<\/h2>\n
Before differential expression analysis, examine your data to determine if there are any confounding issues. Trinity comes with a 'PtR' script that we use to simplify making various charts and plots based on a matrix input file. Run these three commands lines : <\/p>\n
# go to salmon results repertory\ncd \/scratch\/formationX\/ALIGN_AND_ABUNDANCE\/salmon_outdir\n\n#create design file from samples.txt\ncut -f1,2 ..\/samples.txt > design.txt\n\n[orjuela@node25 salmon_outdir]$ more design.txt \nCENPK CENPK_rep1\nCENPK CENPK_rep2\nCENPK CENPK_rep3\nBatch Batch_rep1\nBatch Batch_rep2\nBatch Batch_rep3<\/code><\/pre>\nRun PtR scripts<\/p>\n
\n$path_to_trinity\/Analysis\/DifferentialExpression\/PtR --matrix Trinity_trans.isoform.counts.matrix --samples design.txt --log2 --min_rowSums 10 --compare_replicates\n\n$path_to_trinity\/Analysis\/DifferentialExpression\/PtR --matrix Trinity_trans.isoform.counts.matrix --samples design.txt --log2 --min_rowSums 10 --CPM --sample_cor_matrix\n\n$path_to_trinity\/Analysis\/DifferentialExpression\/PtR --matrix Trinity_trans.isoform.counts.matrix --samples design.txt --log2 --min_rowSums 10 --CPM --center_rows --prin_comp <\/code><\/pre>\nTRANSFERT : Observe plots. Remember, you have to transfert *.pdf files to your home before of transfering into your local machine. <\/p>\n
#transfering plots\ncp *.pdf \/home\/formationX\/\n# from your local machine\nscp formationX@bioinfo-nas.ird.fr:*.pdf .<\/code><\/pre>\n4.2 Identify differentially expressed genes between the two tissues.<\/h2>\n
The tool
run_DE_analysis.pl<\/code> is a PERL script that useBioconductor package edgeR<\/code>. <\/p>\n\n#run DE analysis\n$path_to_trinity\/Analysis\/DifferentialExpression\/run_DE_analysis.pl \\\n--matrix Trinity_trans.isoform.counts.matrix \\\n--method edgeR \\\n--samples_file design.txt \\\n--output edgeR_results<\/code><\/pre>\nThe output files are in the directory edgeR_results. Observe file
Trinity_trans.isoform.counts.matrix.Batch_vs_CENPK.edgeR.DE_results<\/code>. It provides several values for each gene: <\/p>\n1) FDR to indicate whether a gene is differentially expressed or not <\/p>\n
2) logFC is the log2 transformed fold change between the two tissues<\/p>\n
3) logCPM is the log2 transformed normalized read count of average of the samples.<\/p>\n
Usually you have to filter this list of genes\/isoforms to FDR <0.05 or below. To be more conservative, you could also use more stringent FDR cutoff (e.g. <0.001), and only keep genes with high logFC (e.g. <-2 and >2) and\/or high logCPM (e.g. >1). In the edgeR_results directory there is also a \u201cvolcano plot\u201d to visualize the distribution of the DE genes.<\/p>\n
[orjuela@node25 salmon_outdir]$ head edgeR_results\/Trinity_trans.isoform.counts.matrix.Batch_vs_CENPK.edgeR.DE_results\nsampleA sampleB logFC logCPM PValue FDR\nTRINITY_DN332_c0_g1_i10 Batch CENPK -10.359903780639 8.33512668205486 0 0\nTRINITY_DN730_c0_g2_i1 Batch CENPK -6.4740873961699 8.70152404919916 0 0\nTRINITY_DN741_c0_g1_i1 Batch CENPK -6.31231314146213 10.2830859720565 0 0\nTRINITY_DN287_c1_g1_i1 Batch CENPK -6.26650019458591 8.47203996977476 0 0\nTRINITY_DN65_c0_g1_i3 Batch CENPK -4.13189957112457 10.730475730689 0 0\nTRINITY_DN1298_c0_g1_i1 Batch CENPK -3.19109633696326 8.91247526735535 0 0\nTRINITY_DN1253_c0_g1_i1 Batch CENPK -3.03997649547823 8.99502681145441 6.47317750883953e-301 3.75444295512693e-298\nTRINITY_DN51_c0_g1_i1 Batch CENPK -4.59451625218422 10.4521118447878 5.94332889429699e-289 3.01623941385572e-286\nTRINITY_DN708_c0_g1_i1 Batch CENPK -4.1838974791094 7.70916912374135 3.84699757684993e-275 1.73542335133453e-272\n<\/code><\/pre>\n#transfering plots\ncd edgeR_result\ncp *.pdf \/home\/formationX\/\n# from your local machine\nscp formationX@bioinfo-nas.ird.fr:*.pdf .<\/code><\/pre>\n![\"\"](\"https:\/\/southgreenplatform.github.io\/trainings\/\/images\/maplot.png\")
\n![\"\"](\"https:\/\/southgreenplatform.github.io\/trainings\/\/images\/volcano.png\")
<\/p>\n4.3 Clustering analysis<\/h2>\n
Hierarchical clustering and k-means clustering for samples and genes can be done using trinity scripts. <\/p>\n
The clustering will be performed only on differentially expressed genes, with FDR and logFC cutoff defined by -P and -C parameters.<\/p>\n
In this example, we set K=4 for k-means analysis. The genes will be separated into 4 groups based on expression pattern. <\/p>\n
There are two prefiltered files produced:
*DE_results.P1e-3_C2.Batch-UP.subset<\/code> and*DE_results.P1e3_C2.CENPK-UP.subset<\/code>, with differentially expressed genes (FDR cutoff 0.001, logFC cutoff 2 and -2). <\/p>\n#go to edgeR_results\ncd edgeR_results\n# running analyze_diff_expr\n$path_to_trinity\/Analysis\/DifferentialExpression\/analyze_diff_expr.pl \\\n--matrix ..\/Trinity_trans.isoform.TMM.EXPR.matrix --samples ..\/design.txt -P 1e-3 -C 2 \\\n--output cluster_results\n# running define_clusters_by_cutting_tree\n$path_to_trinity\/Analysis\/DifferentialExpression\/define_clusters_by_cutting_tree.pl \\\n-K 4 -R cluster_results.matrix.RData <\/code><\/pre>\n#transfering plots\ncp *.pdf \/home\/formationX\/\n# from your local machine\nscp formationX@bioinfo-nas.ird.fr:*.pdf .<\/code><\/pre>\n![\"\"](\"https:\/\/southgreenplatform.github.io\/trainings\/\/images\/k4.png\")
<\/p>\nBefore finish ...<\/h3>\n
Go to here<\/a> to recover the scripts used in this training.<\/p>\n
\nLicense<\/h3>\n
<\/a><\/p>\n
\nThe resource material is licensed under the Creative Commons Attribution 4.0 International License (here<\/a>).
\n![](\"http:\/\/creativecommons.org.nz\/wp-content\/uploads\/2012\/05\/by-nc-sa1.png\"\/)
\n<\/center>\n<\/div>\n","protected":false},"excerpt":{"rendered":"
- \n
- \n
- Upon successful completion of Trinity, the assembled transcriptome is written to the FASTA file called
- Recover assembly results (only Trinity.fasta) generated by trinity from the shared location to your scratch directory : <\/li>\n<\/ul>\n
- \n
- When the files transfer is finished, verify by listing the content of the current directory and the subdirectory RAWDATA with the
- Copy the exercise files from the shared location to your scratch directory (it is essential that all
- Practice 1:# 1. Check Reads Quality<\/a><\/li>\n