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).

4. Orthology assignment

Given the massive number of genes in the transcriptomes, I primarily used OrthoFinder to assign orthology. In places where orthology remained unclear, I followed up with reciprocal BLASTs to high quality references like A. thaliana and Z. mays on Ensemble Plants and made some gene trees with sequences from Ferrari et al. (2020).

4.1 Downloading the data

I downloaded peptide sequences from a representative set of high quality and of-interest angiosperm genomes for use in OrthoFinder. Most came from Phytozome. Here's a list of species:

Ananas comosus
Amaranthus hypochondriacus
Arabidopsis thaliana
Brachypodium distachyon
Helianthus annuus
Kalanchoe fedtschenkoi
Phalaenopsis equestris
Portulaca amilis (genome proteins)
Portulaca oleracea (TransDecoder longests ORFs)
Sedum album
Sorghum bicolor
Vitis vinifera
Zea mays

(base) [isg4@farnam2 ~]$ cd scratch60/drought-diff-exp/InDepth-Annotation-Pipeline/
(base) [isg4@farnam2 InDepth-Annotation-Pipeline]$ mkdir Reference-proteins
(base) [isg4@farnam2 InDepth-Annotation-Pipeline]$ cd Reference-proteins/
(base) [isg4@c21n06 Reference-proteins]$ curl --cookie jgi_session=/api/sessions/d78dcfb9e60decd737a156a70eb4a378 --output download.20200824.093934.zip -d '{"ids":["591103567ded5e1e49ffd1eb","591103567ded5e1e49ffd1ea","56901a500d878508e3d1fb68","56981cdc0d87851ee9727d1e","564bd7330d878531d71e63b1","564bd7330d878531d71e63ac","569335620d87851ee9726ac8","569335620d87851ee9726ac7","56aef01f0d878559e2870404","5854743e7ded5e78cff8c477","5a3972dc7ded5e35e94f8a0b","5a3972dc7ded5e35e94f8a0c","5e18c789e08d44553ef72280","5e18c789e08d44553ef72286","5e18c782e08d44553ef72207","5e18c782e08d44553ef72212","57bcd2927ded5e0c8712b081","57bcd2947ded5e0c8712b086","5f4058499a211ae42a1a2c48","5f4058499a211ae42a1a2c4a","5bfc629946d1e61e989e5708","5bfc629a46d1e61e989e570a","59cd36857ded5e2f18695ad4","59cd36847ded5e2f18695ad3","582caaed7ded5e2d305b2ad4","582caaec7ded5e2d305b2ad2","582caaf17ded5e2d305b2ae0","582caaf27ded5e2d305b2ae2","5693356b0d87851ee9726afe","5693356b0d87851ee9726af6","5d94dc9ec0d65a87debccfcb","5d94dc9ec0d65a87debccfcc","54ad8ddc0d8785565d4707c5","54ad8ddc0d8785565d4707c2","585473df7ded5e78cff8c43c","585473de7ded5e78cff8c43a","592f5bc37ded5e4e5bbb8e95","592f5bc37ded5e4e5bbb8e97","5cf9223dee2d26b76f024762","5cf9223fee2d26b76f024767","5bfddb7146d1e61e989e789b","5bfddb7046d1e61e989e7899","5c83121646d1e64ae3ba2bec","5c83121546d1e64ae3ba2be8"]}' -H "Content-Type: application/json" https://files.jgi.doe.gov/filedownload/
(base) [isg4@c21n06 Reference-proteins]$ unzip download.20200824.093934.zip
(base) [isg4@c21n06 Reference-proteins]$ mkdir -p Phalaenopsis_equestris/latest_assembly_versions/GCF_001263595.1_ASM126359v1/
(base) [isg4@c21n06 Reference-proteins]$ curl -0 https://ftp.ncbi.nlm.nih.gov/genomes/refseq/plant/Phalaenopsis_equestris/latest_assembly_versions/GCF_001263595.1_ASM126359v1/GCF_001263595.1_ASM126359v1_protein.faa.gz --output Phalaenopsis_equestris/latest_assembly_versions/GCF_001263595.1_ASM126359v1/GCF_001263595.1_ASM126359v1_protein.faa.gz

And I downloaded the [Sedum album] proteins (54539-CDS-prot.fasta) from CoGe.

## YALE FARNAM
(base) [isg4@c21n06 Reference-proteins]$ mkdir Sedum_album

## LOCAL MACHINE
(base) ➜  ~ $ rsync -av ~/Downloads/54539-CDS-prot.fasta isg4@farnam.hpc.yale.edu:/home/isg4/scratch60/drought-diff-exp/InDepth-Annotation-Pipeline/Reference-proteins/Sedum_album

Finally I made a directory of all of the primary proteins in unzipped format.

(base) [isg4@c21n06 Reference-proteins]$ mkdir Primary-proteins
(base) [isg4@c21n06 Reference-proteins]$ find . -name '*.gz' -exec gunzip {} \;
(base) [isg4@c21n06 Reference-proteins]$ mv */*/*/*/*.fa Primary-proteins/
(base) [isg4@c21n06 Reference-proteins]$ mv Phalaenopsis_equestris/latest_assembly_versions/GCF_001263595.1_ASM126359v1/GCF_001263595.1_ASM126359v1_protein.faa Primary-proteins/
(base) [isg4@c21n06 Reference-proteins]$ mv Sedum_album/54539-CDS-prot.fasta Primary-proteins/
(base) [isg4@c21n06 Reference-proteins]$ cp ~/scratch60/Pamilis_funannotate/TRAIN/Portulaca-amilis.v0-FA1.6.0/annotate_misc/genome.proteins.fasta ./Primary-proteins/
(base) [isg4@c21n06 Reference-proteins]$ cp ~/scratch60/drought-diff-exp/de-novo-assembly/P-oleracea-F-no-filter/TransDecoder/Trinity.cluster98.fasta.transdecoder_dir/longest_orfs.pep ./Primary-proteins/

6.2 Running OrthoFinder

(base) [isg4@c21n06 InDepth-Annotation-Pipeline]$ conda activate phylo
(phylo) [isg4@c21n06 InDepth-Annotation-Pipeline]$ mkdir OrthoFinder
(phylo) [isg4@c21n06 InDepth-Annotation-Pipeline]$ cd OrthoFinder
(phylo) [isg4@c21n06 OrthoFinder]$ module load OrthoFinder

Then I ran the following batch file:

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

orthofinder.py -f ../Reference-proteins/Primary-proteins/ -t 20

5. Gene co-expression network analysis

The bulk of this work is in R and Python notebooks. Specifically, WGCNA-P-amilis.ipynb and WGCNA-P-oleracea.ipynb contain the WGCNA analyses that filter genes and libraries, optimize network parameters, and build networks. PA-Network-plotting.ipynb and PO-network-plotting.ipynb contain the subsequent analyses on the subsetted photosynthetic gene networks.

6. Regulatory motif analysis

After identifying different gene modules, we can look back to the genome to see if there are any cis-regulatory elements that may be related to the expression profile recovered. The MEME suite has many, somewhat overlapping tools, which can make it difficult to pick a particular analysis. Our goal is to find if there are shared regulatory elements that may be linked to C4 or CAM expression patterns. In order to test this I first created a de novo set of Portulaca motifs and downloaded an existing curated database of plant motifs. The curated motifs come from JASPAR (Khan et al. 2018). Then I removed redundant motifs between the de novo and JASPAR database, and tested for motif enrichment in multiple modules of genes from WGCNA.

6.1 de novo motif identification

6.1.1 Extracting genomic regions

In order to identify putative regulatory motifs from the P. amilis genome, I first had to extract the various regions from the genomic scaffolds. In the notebook Regulatory-element-extractor.ipynb I used gffutils to parse gff and tbl files of genomic features and extract regions from the genomic scaffolds. These include

26,969 5'UTRs
208,192 introns
25,336 3'UTRs
57,188 2kb upstream windows

These files have MEME-recognizable fasta headers, for example

>scaffold_1:2298-2735 ID=FUN_000001-T1.exon4_1;Strand=(-);Source=gffutils_derived

This header can be parsed by MEME programs so that the resulting gff files have positions that are relative to the scaffolds instead of relative to the feature they are derived from.

8.2.2 Extracting relevant sequences

Now that I've gotten these non-coding regions properly formatted, I need a way of extracting sequences. As Bailey (2007) describes, motif discovery is most successful when you have well chosen inputs. He reccomends:

Include as many sequences as possible that contain the motifs.

Keep the sequences as short as possible.

Remove sequences that are unlikely to contain any motifs.

I won't be altering the lengths of the sequences here, but we can leverage our co-expression results to search for motifs withinn co-expression modules. To do this I wrote a short script that takes a file containing the list of gene models (one per line) and a fasta to extract from. Gene-model-extractor.py pulls out any records that match the gene models. I'll begin by getting gene lists for each WGCNA module for both water and drought.

I've created a list of all of the photosynthetic genes (phang-genemodels.txt) and placed it in a new Motif-analysis directory on my local machine. Here's an example of how to extract all of the 5'UTRs for the photosynthetic gene models.

## LOCAL MACHINE
(phylo) ➜  Drought_Experiment_3 git:(master) ✗ $ mkdir -p Motif-analysis
(phylo) ➜  Drought_Experiment_3 git:(master) ✗ $ python ./Scripts/Gene-model-extractor.py --models ./Motif-analysis/phang-genemodels.txt --fasta ~/Dropbox/GitHub_repos/finalproject/scratch/Pamilis_funannotate/TRAIN/Portulaca-amilis.v0-FA1.6.0/update_results/Portulaca_amilis.five_prime_UTR.FIMO.fasta --output ./Motif-analysis/Gene-sets/P-amilis.phangs.5UTR.FIMO.fasta

   _____ __  __ ______
  / ____|  \/  |  ____|
 | |  __| \  / | |__
 | | |_ | |\/| |  __|
 | |__| | |  | | |____
  \_____|_|  |_|______|

Gene Model Extractor

Created output directory /Users/iangilman/Dropbox/Drought-experiments/Drought_Experiment_3/Motif-analysis/Gene-sets
Model labels path: Motif-analysis/phang-genemodels.txt
Master fasta path: /Users/iangilman/Dropbox/GitHub_repos/finalproject/scratch/Pamilis_funannotate/TRAIN/Portulaca-amilis.v0-FA1.6.0/update_results/Portulaca_amilis.five_prime_UTR.FIMO.fasta
Output path: Motif-analysis/Gene-sets/P-amilis.phangs.5UTR.FIMO.fasta
Portulaca_amilis.five_prime_UTR.FIMO.fasta contains 26969 records
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 453/453 [00:04<00:00, 97.91it/s]
[04 Jan 2021 15:02:45] Time elapsed: 0:00:05.025076

Now that the script is up and running, I'll extract each well-watered and drought set of genes. I first used a little CLI iPython to generate these lists and deposited them in Motif-analysis/Gene-sets, and then ran bash (I actually run zsh on my own computer) loops to extract the relvant features. For example,

(phylo) ➜  Gene-sets git:(master) ✗ $ for f in Drought-modules/*.txt;
do;
NAME=$f:t:r;
echo ./Drought-modules/${NAME}.fasta;
python ../../Scripts/Gene-model-extractor.py --models $f --fasta ~/Dropbox/GitHub_repos/finalproject/scratch/Pamilis_funannotate/TRAIN/Portulaca-amilis.v0-FA1.6.0/update_results/Portulaca_amilis.five_prime_UTR.FIMO.fasta --output ./Drought-modules/${NAME}.5UTR.fasta;
done

I repeated this with 3'UTR, 2kb upstream promoter, and intronic regions and then sent the results to the cluster.

(phylo) ➜  Motif-analysis git:(master) ✗ $ rsync -av Gene-sets isg4@farnam.hpc.yale.edu:/home/isg4/scratch60/drought-diff-exp/Motif-analysis

I forgot I had to do this for the intersection ppc-1E1a' and ppc-1E1c drought and well watered modules. I pulled those lists from the Network-plotting notebook and then extracted the respective regions.

6.1 Motif databases

The plant specific JASPAR CORE database can be downloaded in MEME format here.

6.2 Genomic region extraction

I first wrote a Python script that takes a set of genomic scaffolds and an associated gff file and extracts the possible regulatory regions. I considered a 2000bp upstream region (putative promoter), intronic regions, 5' UTRs, and 3'UTRs. After generating these sequence files I used the MEME suite to detect enriched elements of interest.