Home

In this Wiki I'll walk through the transcription-related analyses of Portulaca amilis and P. oleracea that formed the basis for our recent paper

Gilman IS, Moreno-Villena J, Lewis ZR, Goolsby EW, Edwards EJ. submitted. Gene co-expression networks reveal orthology among multiple photosynthetic pathways in Portulaca

that investigated the two photosynthetic systems of Portulaca (C4 and CAM). I'll try to include everything necessary to reproduce our results and figures.

The two main sources of information we are using to investigate Portulaca's C4+CAM system are 1) the P. amilis genome and 2) transcripts from a CAM induction experiment. All raw RNAseq data can be found on the NCBI's SRA under BioProject PRJNA732408, and the P. amilis genome is hosted on Phtyozome. These data will be used separately to understand the transcript abundance of genes and their genomic context in the P. amilis genome. The P. amilis transcriptome was also used during the genome annotation process, but in this wiki I'll stick to

P. oleracea de novo transcriptome assembly
transcript differential abundance analyses (aka 'differential expression')
gene co-expression network analysis
cis-element analysis

I've found this review by Babarinde et al. (2019) to be very helpful in understanding the mechanics and design of RNAseq projects. References within that review, as well as this paper by Wilhelmsson et al. (2019) highlight that there are trade offs between different methods, and that employing a few different methods or a consensus of results is best if possible. During the evolution of this paper I worked through many transcriptome assembly processes with varying filtering and trimming stringencies, multiple read quantification tools (RSEM, STAR, Kallisto), multiple clustering techniques (maSigPro, WGCNA, and custom methods), and multiple quantification programs (EdgeR, DESeq2, and sleuth), but I'll mostly show the final analyses here.

0. Setting up

I'll be continuing with the super_funannotate conda environment that I used during genome annotation, but note that I change environments often. You can always see what environment I'm using by looking at my command line prefix, for example, below I change from the base to super_funannotate environment.

(base) [isg4@c19n09 scratch60]$ conda activate super_funannotate
(super_funannotate) [isg4@c19n09 scratch60]$ mkdir drought-diff-exp
(super_funannotate) [isg4@c19n09 scratch60]$ cd drought-diff-exp

I'd already downloaded all of sequence data from the Yale Center for Genome Analysis servers into our shared storage space, but I made a working copy of those directories in a working data directory.

(super_funannotate) [isg4@c19n09 scratch60]$ mkdir data
(super_funannotate) [isg4@c19n09 drought-diff-exp]$ cp -r path/to/my/data/Portulaca/Drought_Exp-2019/P_amilis/* ./data/
(super_funannotate) [isg4@c19n09 drought-diff-exp]$ cp -r path/to/my/data/Portulaca/Drought_Exp-2019/P_oleracea-Florida/ ./data/

1. Quality control and read trimming

The first thing I did with the raw reads was look at a few quality metrics with FastQC and filtered and trimmed reads with PRINSEQ. I tried using Trimmomatic but found it wasn't removing PolyA/T tails properly (more below).

1.1 Quality control with FastQC

FastQC analyses metrics like per base/sequence quality, GC content, N content, sequence length distribution, duplication levels, and adapter content. It nicely summarizes these in html files with visuals to discern between low, moderate, and high quality data.

(super_funannotate) [isg4@c14n03 drought-diff-exp]$ mkdir -p quality_control/FastQC
(super_funannotate) [isg4@c14n03 drought-diff-exp]$ cd quality_control/

Then, in an array-formatted batch file, I ran my FastQC command:

fastqc -t 12 --outdir ./FastQC/ ../data/*/*.fastq.gz

I found that my forward (_R1_) reads frequently had poly-A as an over-represented sequence.

1.2 Quality control and trimming with PRINSEQ

Here's an example of a PRINSEQ run on one of the P. oleracea files. Note that PRINSEQ cannot take zipped files as input, so I decompressed them first.

(super_funannotate) [isg4@c14n09 drought-diff-exp]$ makdir -p quality_control/PRINSEQ/P-oleracea-F_001_186_007
(super_funannotate) [isg4@c14n09 drought-diff-exp]$ cd quality_control/PRINSEQ
(super_funannotate) [isg4@c14n09 PRINSEQ]$ for f in ../../data/P-oleracea-F_001_186_007/*.gz; do STEM=$(basename "${f}" .gz); echo $STEM; gunzip -c "${f}" > ./P-oleracea-F_001_186_007/"${STEM}";done
(super_funannotate) [isg4@c14n09 PRINSEQ]$ cd P-oleracea-F_001_186_007/
(super_funannotate) [isg4@c14n09 P-oleracea-F_001_186_007]$ prinseq-lite.pl -verbose -fastq ISG_001_186_007_S17_L001_R1_001.fastq -fastq2 ISG_001_186_007_S17_L001_R2_001.fastq -out_format 3 -min_qual_mean 25 -trim_tail_left 12 -trim_tail_right 12 -trim_qual_left 20 -trim_qual_right 20 -ns_max_p 1  -stats_all

Here is an explanation the flags:

min_qual_mean 25 : filter sequences with mean Phred < 25
trim_tail_left/right 12 : trim poly-A/T tails longer than 12bp
trim_qual_left/right 20 : trim tails with Phred < 20
ns_max_p 1 : filter out sequences with 1% or more Ns

This run ended up removing about 250,000 reads (~1%) from this library because they were singletons, had low quality scores, or were too short after trimming. PRINSEQ did a good job removing the poly-A/T tails (found in ~34,000 reads). The remaining libraries were processed with the following batch file, which unzipped .fastq.gz files, ran PRINSEQ, and then produced the related graphs based on the output

#!/bin/bash
#SBATCH --job-name=prinseq
#SBATCH --time=12:00:00
#SBATCH -c 5
#SBATCH -p general
#SBATCH --mem=24G
#SBATCH --array=1-92
#SBATCH --ntasks-per-node=1
#SBATCH --mail-type=ALL
#SBATCH --output=prinseq-%A-%a.out
#SBATCH --error=prinseq-%A-%a.err

module load PRINSEQ

sleep $((SLURM_ARRAY_TASK_ID*60))

SAMPLES=(
P-oleracea-F_005_138_055
P-oleracea-F_011_161_032
P-oleracea-F_013_137_056
...
)

SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID-1]}
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Begin $SAMPLE"

mkdir -p ./quality_control/PRINSEQ/${SAMPLE}
cd ./quality_control/PRINSEQ/${SAMPLE}

dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Unzipping"

for f in ~/scratch60/drought-diff-exp/data/${SAMPLE}/*.gz
        do STEM=$(basename "${f}" .gz)
	echo $STEM
        gunzip -c "${f}" > ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}/"${STEM}"
done

cd ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}

dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Begin PRINSEQ"

prinseq-lite.pl -verbose -fastq *_R1_001.fastq -fastq2 *_R2_001.fastq -out_good ${SAMPLE}.G -out_bad null -log ${SAMPLE}.log -graph_data ${SAMPLE}.gd -out_format 3 -min_qual_mean 25 -trim_tail_left 12 -trim_tail_right 12 -trim_qual_left 20 -trim_qual_right 20 -ns_max_p 1

dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Finished PRINSEQ"

for g in ./*.gd
        do STEM2=$(basename "${g}" .gd)
        dt=$(date '+%d/%m/%Y %H:%M:%S')
	echo [$dt] "Begin graphing $STEM2"
	prinseq-graphs.pl -log -html_all -png_all -i ${STEM2}.gd -o ${STEM2}.gd
done

dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Finished graphing"

2. Transcriptome assembly

The P. amilis transcriptome assembly took place during the annotation of the P. amilis genome (see funannotate train section), but ultimately the P. amilis reads were mapped to the genome, and the transcriptome was primarily being used for annotation purposes. For P. oleracea, I assembled a de novo transcriptome using trinity. I began by creating concatenated read files of all of the P. oleracea samples. This was probably overkill and we could've probably just skimmed ~10M reads off the top of each file (instead of all ~20-25M), but greater read counts translates into greater k-mer counts, and in turn to an increased ability to call rare isoforms.

2.1 Trinity de novo assembly

I tried a few iterations of Trinity assemblies but ultimately ended up with a pretty vanilla methodology. I initially tried filtering rRNA using SortMeRNA, which removed ~5-15% of reads per sample, but resulted in an assembly with log contig sizes and terrible BUSCOs scores. I also tried (for posterity) running a P. amilis genome-guided run of Trinity, but that yielded even smaller contigs and still worse BUSCO scores. In the end I went with a simple assembly of the reads output from PRINSEQ.

I gave each set of reads to Trinity in line, and ran the latest version (2.11.0) in a singularity container within a batch file.

singularity exec -e trinityrnaseq.v2.11.0.simg Trinity --max_memory 500G --seqType fq --output ./P-oleracea-F-no-filter/trinity --CPU 20 --trimmomatic --left \
/home/isg4/scratch60/drought-diff-exp/quality_control/PRINSEQ/P-oleracea-F_001_186_007/P-oleracea-F_001_186_007.G_1.fastq,/home/isg4/scratch60/drought-diff-exp/quality_control/PRINSEQ/P-oleracea-F_005_138_055/P-oleracea-F_005_138_055.G_1.fastq,... \
--right \
/home/isg4/scratch60/drought-diff-exp/quality_control/PRINSEQ/P-oleracea-F_001_186_007/P-oleracea-F_001_186_007.G_2.fastq,/home/isg4/scratch60/drought-diff-exp/quality_control/PRINSEQ/P-oleracea-F_005_138_055/P-oleracea-F_005_138_055.G_2.fastq,...

2.2 Redundancy reduction

I used CD-HIT to cluster contigs with 98% similarity or greater.

(phylo) [isg4@p02r09n20 P-oleracea-F-no-filter]$ mkdir CD-HIT/
(phylo) [isg4@p02r09n20 P-oleracea-F-no-filter]$ cd CD-HIT/
(phylo) [isg4@p02r09n20 CD-HIT]$ cd-hit-est -o Trinity.cluster98 -c 0.98 -i ../trinity/Trinity.fasta -p 1 -d 0 -b 3 -T 0

This reduced the number of transcripts from 444,307 to 413,658. So not a huge reduction.

################################
Total trinity 'genes':  230086
Total trinity transcripts:      413658
Percent GC: 40.58

########################################
Stats based on ALL transcript contigs:
########################################

        Contig N10: 3757
        Contig N20: 2817
        Contig N30: 2251
        Contig N40: 1820
        Contig N50: 1446

        Median contig length: 520
        Average contig: 882.41
        Total assembled bases: 365016900


#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

        Contig N10: 3507
        Contig N20: 2521
        Contig N30: 1883
        Contig N40: 1385
        Contig N50: 1019

        Median contig length: 394
        Average contig: 684.43
        Total assembled bases: 157478620

So we removed isoforms almost exclusively using CD-HIT, as the the number of "genes" barely changed. The length statistics slightly decreased, meaning that we removed some longer transcripts.

2.3 Transcriptome quality assessment

The Trinity quality assessment page has a great list of quality control suggestions, some of which I show below, and I already presented some of the Trinity length statistics, above.

2.3.1 Assessing read content

I assessed the number of reads mapping to the transcriptome using kallisto and following the steps outlined here.

(funannotate) [isg4@p02r09n20 drought-diff-exp]$ cd Kallisto/
(funannotate) [isg4@p02r09n20 drought-diff-exp]$ mkdir -p P-oleracea/no-filter-cluster98

In a batch file I ran the kallisto index command

kallisto index -i P-oleracea-F-no-filter-cluster98-transcript.idx ~/scratch60/drought-diff-exp/de-novo-assembly/P-oleracea-F-no-filter/CD-HIT/Trinity.cluster98.fasta

and then read quantification in an array job

#!/bin/bash
#SBATCH --job-name=kallisto-array
#SBATCH --time=6:00:00
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=24G
#SBATCH --array=1-44
#SBATCH --ntasks-per-node=1
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=%x-%A-%a.out
#SBATCH --error=%x-%A-%a.err

sleep $((SLURM_ARRAY_TASK_ID*60))

SAMPLES=(
P-oleracea-F_001_186_007
P-oleracea-F_005_138_055
...
P-oleracea-F_142_127_066
)

SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID-1]}
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Begin $SAMPLE"

kallisto quant -i P-oleracea-F-no-filter-cluster98-transcript.idx -o $SAMPLE -b 100 --threads 20 /home/isg4/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}/${SAMPLE}.G_1.fastq /home/isg4/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}/${SAMPLE}.G_2.fastq

dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Finished $SAMPLE"

Checking the resulting json files, 63.7+/-2.4% of reads mapped back to the transcriptome per sample. This isn't horrible but it isn't great either. Ideally we would have been around 80% or better.

2.3.2 BUSCO

I also ran BUSCO v5 through a singularity container for both embryophyta and eukaryota.

(base) [isg4@p02r09n20 P-oleracea-F-no-filter]$ cd BUSCO/
(base) [isg4@p02r09n20 BUSCO]$ mkdir P-oleracea-F-no-filter-cluster98
(base) [isg4@p02r09n20 BUSCO]$ cd P-oleracea-F-no-filter-cluster98/
(base) [isg4@p02r09n20 P-oleracea-F-no-filter-cluster98]$ singularity exec -e docker://ezlabgva/busco:v5.beta.1_cv1 busco -i ../../CD-HIT/Trinity.cluster98.fasta -o cluster98-embryophyta -m transcriptome -l embryophyta -c 20 --out_path ./embryophyta
(base) [isg4@p02r09n20 P-oleracea-F-no-filter-cluster98]$ singularity exec -e docker://ezlabgva/busco:v5.beta.1_cv1 busco -i ../../CD-HIT/Trinity.cluster98.fasta -o cluster98-eukaryota -m transcriptome -l eukaryota -c 20 --out_path ./eukaryota

	Embryophyta (1,614)	Eukaryota (255)
Complete	1,527 (94.7%)	252 (98.9%)
Single copy	264 (16.4%)	43 (16.9%)
Duplicated	1,263 (78.3%)	209 (82.0%)
Fragmented	52 (3.2%)	0 (0.0%)
Missing	35 (2.1%)	3 (1.1%)

3. Transcript abundance

3.1 Read quantification

Read quantification was already done for P. oleracea as part of 2.3 Assessing read content. The P. amilis analyses followed the same structure, but the index file came from the update step of funannotate.

(super_funannotate) [isg4@c15n07 drought-diff-exp]$ mkdir -p Kallisto/Portulaca-amilis
(super_funannotate) [isg4@c15n07 drought-diff-exp]$ cd Kallisto/Portulaca-amilis
(super_funannotate) [isg4@c15n07 Portulaca-amilis]$ kallisto index -i ./P-amilis-update-mRNA.index ~/scratch60/Pamilis_funannotate/TRAIN/Portulaca-amilis.v0-FA1.6.0/update_results/Portulaca_amilis.mrna-transcripts.fa > P-amilis-update-mRNA.index.log 2>&1

Then a similar batch script was used to analyze all P. amilis samples from PRINSEQ.

#!/bin/bash
#SBATCH --job-name=kallisto
#SBATCH --time=6:00:00
#SBATCH -c 20
#SBATCH -p general
#SBATCH --mem=24G
#SBATCH --array=1-48
#SBATCH --ntasks-per-node=1
#SBATCH --mail-type=ALL
#SBATCH --mail-user=ian.gilman@yale.edu
#SBATCH --output=kallisto-%A-%a.out
#SBATCH --error=kallisto-%A-%a.err

sleep $((SLURM_ARRAY_TASK_ID*60))

SAMPLES=(
P-amilis_002_174_019
P-amilis_004_150_043
...
P-amilis_096_295_282
)

SAMPLE=${SAMPLES[$SLURM_ARRAY_TASK_ID-1]}
dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Begin $SAMPLE"

kallisto quant -i P-amilis-update-mRNA.index -o $SAMPLE -b 100 --threads 20 --genomebam --gtf ~/scratch60/drought-diff-exp/RSEM/Portulaca-amilis/Portulaca-amilis.gtf ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}/${SAMPLE}.G_1.fastq ~/scratch60/drought-diff-exp/quality_control/PRINSEQ/${SAMPLE}/${SAMPLE}.G_2.fastq

dt=$(date '+%d/%m/%Y %H:%M:%S')
echo [$dt] "Finished $SAMPLE"

3.2 Differential abundance analysis with sleuth

The R notebook with the analyses can be found here, so I'll only briefly describe the major points, and I followed this vignette when creating the analyses. Using kallisto + sleuth allows us to account for uncertainty from the very beginning of the analysis (pseuoaligning reads) and better partition biological and technical variance. These methods also maintain the compatibility between transcript- and gene-level counts through the use of a transcript-gene map and p-value aggregation.

First, in order to associate genes with transcripts, I'll grab all of the CDS from the annotation.

(super_funannotate) [isg4@c14n08 update_results]$ grep ">" Portulaca_amilis.cds-transcripts.fa > transcript-gene-map.txt
(super_funannotate) [isg4@c14n08 update_results]$ sed -i -e 's/>//g' transcript-gene-map.txt

## LOCAL MACHINE
(base) ➜  farnam-drought-diff-exp git:(master) ✗ $ rsync -av isg4@farnam.hpc.yale.edu:/Portulaca-amilis.v0-FA1.6.0/update_results/transcript-gene-map.txt .

From there, the P. amilis and P. oleracea analyses proceeded similarly. The heart of this analysis lies in the comparison between the fits of the full model (counts ~ time_point + treatment) and the reduced model (counts ~ time_point).

## Time only
so = sleuth_fit(obj = so, formula = ~time_point, fit_name = 'reduced')
## Full model
so = sleuth_fit(obj = so, formula = ~time_point + treatment, fit_name = 'full')
## Compare
so = sleuth_lrt(obj = so, null_model = 'reduced', alt_model = 'full')

This likelihood ratio test identifies those genes who patterns of expression are better explained when accounting for the treatment (well-watered vs. drought-stressed).