Phaeocystis globosa colonial gene expression

<p>Data and analysis for the&nbsp;paper:&nbsp;</p> <p><strong>Differential gene expression supports a resource-intensive, defensive role for colony production in the bloom-forming haptophyte, <em>Phaeocystis globosa</em></strong></p> <p>by:&nbsp;Margaret Mars Brisbin and Satoshi Mitarai</p> <p>The <em>Phaeocystis&nbsp;globosa</em> CCMP1528 transcriptome used in the study (phaeocystisglobosa_euk_seqs.fasta&nbsp;or pg_euk_seqs_altnames.fasta) was assembled with trimmed sequencing reads from 8 biological replicates (4 colonial replicates and 4 solitary replicates) with the Trinity software (v2.3.2).</p> <p>Raw sequencing reads are available from the NCBI SRA with accession numbers:&nbsp;SRR7811979&ndash;SRR7811986.</p> <p>Before assembling the transcriptome, reads were quality filtered and trimmed with the Trimmomatic software (v3.36) using the command:</p> <pre><code>java -jar $TRIM/trimmomatic-0.36.jar PE -phred33 $DATA2/S${SLURM_ARRAY_TASK_ID}_S*_R1_001.fastq.gz \ $DATA2/S${SLURM_ARRAY_TASK_ID}_S*_R2_001.fastq.gz \ $OUT/S${SLURM_ARRAY_TASK_ID}_1_paired.fq $OUT/S${SLURM_ARRAY_TASK_ID}_1_unpaired.fq \ $OUT/S${SLURM_ARRAY_TASK_ID}_2_paired.fq $OUT/S${SLURM_ARRAY_TASK_ID}_2_unpaired.fq \ ILLUMINACLIP:$TRIM/adapters/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36</code></pre> <p>Trimmed reads were mapped&nbsp;to the ERCC reference sequences for Mix1 and mapped reads were filtered using the following commands from bowtie2 (v2.2.6), samtools, and bedtools:&nbsp;</p> <pre><code>bowtie2 -t -x $REF \ -1 $DATA/S${SLURM_ARRAY_TASK_ID}_1_paired.fq \ -2 $DATA/S${SLURM_ARRAY_TASK_ID}_2_paired.fq \ -S $OUT/S${SLURM_ARRAY_TASK_ID}_ercc.sam samtools view -bS $DATA/S${SLURM_ARRAY_TASK_ID}_ercc.sam &gt;$DATA/S${SLURM_ARRAY_TASK_ID}.bam samtools sort $DATA/S${SLURM_ARRAY_TASK_ID}.bam $DATA/S${SLURM_ARRAY_TASK_ID}_sorted samtools view -b -f 13 S${SLURM_ARRAY_TASK_ID}_sorted.bam &gt; S${SLURM_ARRAY_TASK_ID}_unmapped.bam samtools sort -n $DATA/S${SLURM_ARRAY_TASK_ID}_unmapped.bam $DATA/S${SLURM_ARRAY_TASK_ID}.qsort bedtools bamtofastq -i $DATA/S${SLURM_ARRAY_TASK_ID}.qsort.bam -fq $DATA/S${SLURM_ARRAY_TASK_ID}_1_paired.fq -fq2 $DATA/S${SLURM_ARRAY_TASK_ID}_2_paired.fq</code></pre> <p>The resulting Trimmed reads&nbsp;without ERCC sequences were used to make the transcriptome assembly:&nbsp;</p> <pre><code>Trinity --seqType fq --max_memory 475G \ --left $DATA2/C1_1_paired.fq,$DATA2/C2_1_paired.fq,$DATA2/C3_1_paired.fq,$DATA2/C4_1_paired.fq,$DATA2/S1_1_paired.fq,$DATA2/S2_1_paired.fq,$DATA2/S3_1_paired.fq,$DATA2/S4_1_paired.fq \ --right $DATA2/C1_2_paired.fq,$DATA2/C2_1_paired.fq,$DATA2/C3_2_paired.fq,$DATA2/C4_2_paired.fq,$DATA2/S1_2_paired.fq,$DATA2/S2_2_paired.fq,$DATA2/S3_2_paired.fq,$DATA2/S4_2_paired.fq \ --CPU 12</code></pre> <p>The Trinity assembly was dereplicated with CD-HIT-EST (v2016-0304) at 95% :&nbsp;</p> <pre><code>cd-hit-est -i $DATA/Trinity.fasta -o Trinity_Pg_clustered_95 -c 0.95 -n 8 -p 1 -g 1 -M 200000 -T 8 -d 40</code></pre> <p>The Trinity assembly was filtered to remove bacterial contamination&nbsp;by first running a blastn(v2.6.0+) against the nr/nt NCBI database:</p> <pre><code>blastn -query $DATA/Trinity_Pg_clustered_95.fasta -task blastn -db $REF -num_threads 12 -max_target_seqs 1 -outfmt 5 &gt; TrinityBlast.xml</code></pre> <p>and then&nbsp;removing bacterial reads&nbsp;with custom python scripts included here: TrinityBlastXML.ipynb and&nbsp;FIlterTrinityEukNotEuk.ipynb&nbsp;</p> <p>RSEM (v1.2.22) was run with the final transcriptome assembly (phaeocystisglobosa_euk_seqs.fasta&nbsp;or pg_euk_seqs_altnames.fasta):&nbsp;</p> <pre><code>rsem-calculate-expression --bowtie2 --paired-end \ $DATA/C${SLURM_ARRAY_TASK_ID}_1_paired.fq \ $DATA/C${SLURM_ARRAY_TASK_ID}_2_paired.fq \ $REF/rsemref_longISO/pg_euks_RSEMref \ $REF/rsemout_longISO/C${SLURM_ARRAY_TASK_ID} rsem-calculate-expression --bowtie2 --paired-end \ $DATA/S${SLURM_ARRAY_TASK_ID}_1_paired.fq \ $DATA/S${SLURM_ARRAY_TASK_ID}_2_paired.fq \ $REF/rsemref_longISO/pg_euks_RSEMref \ $REF/rsemout_longISO/S${SLURM_ARRAY_TASK_ID} </code></pre> <p>The resulting data files are: C*.genes.results and S*.genes.results which were used with DESeq2 in the R environment to analyze different gene expression. The code for these analyses is available in&nbsp;html and R markdown (PhaeoColSol_DE.html,&nbsp;PhaeoColSol_DE.Rmd).&nbsp;</p> <p>The transcriptome assembly was annotated with the Dammit software (v1.0rc2), which wraps Transdecoder, HMMER, and BUSCO, and by submitting the translated amino acid sequences to GhostKOALA.&nbsp;</p> <p>The raw pfam Dammit annotation results are included: pg_euk_seqs.fasta.x.pfam.gff3. These results were parsed with the script: Pfam_gffParsing.ipynb. The resulting file,&nbsp;pfam_parsed_annotation.csv,&nbsp;is used in the script&nbsp;PhaeoColSol_DE.Rmd with pfam2go4R.txt for GO enrichment analysis. The script shinycolsol.Rmd&nbsp;creates an interactive plot of GO enrichment results.&nbsp;</p> <p>The GhostKOALA results are&nbsp;user_ko.csv, and are&nbsp;used in the script&nbsp;PhaeoColSol_DE.Rmd for KEGG pathway enrichment analysis.&nbsp;</p>