genefish 23:40 on April 30

Shelly’s Notebook: Apr 24-30, 2020 Sea lice methylome

Sample details

2 lice samples from 20190523
- 2 adult females
Samples from 20191118 (in fridge in 209)
- 2 pools of adult sea lice from different populations
- from two different salmon farms
- different salinities? temperatures?
*** need more info about both of sample groups; waiting to hear back from Cristian’s lab

Methylation data

re-trimmed all samples with bismark recommended trimming parameters

script here: https://gannet.fish.washington.edu/metacarcinus/mox_jobs/20200427_TG_SalmoCalig.sh
slurm file here: https://gannet.fish.washington.edu/metacarcinus/Salmo_Calig/analyses/20200427/slurm-2535023.out
trimmed reads here: https://gannet.fish.washington.edu/metacarcinus/Salmo_Calig/analyses/20200427/TG_PE_FASTQS/

aligned trimmed reads to genome

script here: https://gannet.fish.washington.edu/metacarcinus/mox_jobs/20200427_Align_Calig.sh
slurm file here: https://gannet.fish.washington.edu/metacarcinus/Salmo_Calig/analyses/20200427/TG_PE_Calig_Aligned/slurm-2537134.out
bismark output here: https://gannet.fish.washington.edu/metacarcinus/Salmo_Calig/analyses/20200427/TG_PE_Calig_Aligned/
alignment summary:
- generated the summary table below with this R markdown using the multiqc results:

desc. stat	Sealice F1 S20	Sealice F2 S22
total reads before trim	226181237	163876559
perc reads trim removed	0.56	0.43
total reads after trim	224907116	163179589
uniq aligned reads	62174415	34864893
perc uniq aligned	27.64	21.37
ambig reads	49967995	31414190
perc ambig aligned	22.22	19.25
perc no align	50.14	59.38
dedup reads	39280061	26425601
dedup reads percent	63.18	75.79
dup reads	22894354	8439292
dup reads percent	36.82	24.21
percent cpg meth	1.1	1.0
percent chg meth	1.0	0.8
percent chh meth	1.4	1.5

qualimap summary of alignments:
- bamqc script: 20200429_sealice_qualimap.sh
- log file: 20200429_sealice_qualimap.log
- output:
- multibamqc script: 20200429_sealice_qualimap2.sh
- log file: 20200429_sealice_qualimap2.log
- output: https://gannet.fish.washington.edu/metacarcinus/Salmo_Calig/analyses/20200429/
  - report: [https://gannet.fish.washington.edu/metacarcinus/Salmo_Calig/analyses/20200429/multisampleBamQcReport.html](https://gannet.fish.washington.edu/metacarcinus/Salmo_Calig/analyses/20200429/multisampleBamQcReport.html)
  - Coverage:
    - horizontal:
      - For F1_S20: ~85% of the genome is covered at 1x, ~65% is covered at 5x
      - F1_S20:
      - For F2_S22: ~75% of the genome is covered at 1x, ~45% is covered at 5x
        F2 got ~168M reads vs F1 which got 226 reads
      - F2_S22:
      - here’s an example of data (human) that is a little more saturated (has more genome coverage at higher depth; more than 90% of the genome is covered at 20x read depth)
    - vertical coverage (depth) of scaffold bases:
      - F1: 16.4619x +/- 61.5992 (s.d.)
      - F2: 10.698x +/- 39.3661 (s.d.)
    - ran bedtools genomecov for comparison and got the same horizontal coverage info
      - jupyter notebook here: https://github.com/shellytrigg/Salmon_sealice/blob/master/jupyter/SealiceGenomeCoverage.ipynb

prepared merged CpG 5x cov files

categorized CpGs with 5x cov:

jupyter notebook: https://github.com/shellytrigg/Salmon_sealice/blob/master/jupyter/20200429_Sealice_mCpG_analysis.ipynb
unmethylated CpGs (< 10% mCpG):
- Sealice_F1_S20.dedup.cov.5x.unmeth.bed
- Sealice_F2_S22.dedup.cov.5x.unmeth.bed
sprasely methylated CpGs (10-50%):
- Sealice_F1_S20.dedup.cov.5x.sparmeth.bed
- Sealice_F2_S22.dedup.cov.5x.sparmeth.bed
methylated CpGs (> 50%):
- Sealice_F1_S20.dedup.cov.5x.meth.bed
- Sealice_F2_S22.dedup.cov.5x.meth.bed
merged unmethylated CpGs (< 10% mCpG):
- Sealice_F1_S20.dedup.cov.5x.merge.unmeth.bed
- Sealice_F2_S22.dedup.cov.5x.merge.unmeth.bed
merged sprasely methylated CpGs (10-50%):
- Sealice_F1_S20.dedup.cov.5x.merge.sparmeth.bed
- Sealice_F2_S22.dedup.cov.5x.merge.sparmeth.bed
merged methylated CpGs (> 50%):
- Sealice_F1_S20.dedup.cov.5x.merge.meth.bed
- Sealice_F2_S22.dedup.cov.5x.merge.meth.bed

5x CpG summary tables:

Sample|Methylated CpG (>= 50%)|Sparsely methylated CpG (10% – 50%)|Unmethylated CpG (< 10%) :—–:|:—–:|:—–:|:—–: F1|2335|391515|8890795 F2|1864|2232274|6108325

CpG category|F1:F2 CpG overlap|Uniq F1 CpGs|Uniq F2 CpGs|frac F1 mCpG overlapping|frac F2 mCpG overlapping :—–:|:—–:|:—–:|:—–:|:—–:|:—–: Methylated CpG(>= 50%) | 545 | 1790 | 1319 |23.34 | 29.24 Sparsely methylated CpG(10% – 50%) | 19838 | 371677 |212436| 5.07 | 8.54 Unmethylated CpG(< 10%) |5431008 |3459787| 677317 |61.09| 88.91

5x merged CpG summary tables:

Sample	Methylated CpG (>= 50%)	Sparsely methylated CpG (10% – 50%)	Unmethylated CpG (< 10%)
F1	1342	233941	6831337
F2	1102	165423	5130433

CpG category	F1:F2 CpG overlap	Uniq F1 CpGs	Uniq F2 CpGs	frac F1 mCpG overlapping	frac F2 mCpG overlapping
Methylated CpG(>= 50%)	314	1028	788	23.40	28.49
Sparsely methylated CpG(10% – 50%)	13570	220371	151853	5.80	8.20
Unmethylated CpG(< 10%)	4824175	2007162	306258	70.62	94.03

IGV session

https://github.com/shellytrigg/Salmon_sealice/blob/master/analyses/20200429_CpG_analysis/20200429_CpG_analysis.xml

Example of highly methylated CpG overlapping between both samples Screen%20Shot%202020-04-29%20at%202.54.54%20PM.png

zoomed in view Screen%20Shot%202020-04-29%20at%202.55.57%20PM.png

Genomic Feature analysis

used the same jupyter notebook: https://github.com/shellytrigg/Salmon_sealice/blob/master/jupyter/20200429_Sealice_mCpG_analysis.ipynb
downloaded gff from here: https://figshare.com/articles/Chromosome-scale_genome_assembly_of_the_sea_louse_Caligus_rogercresseyi/11780658
- file is here: https://gannet.fish.washington.edu/metacarcinus/Salmo_Calig/GENOMES/Caligus/Caligus-rogercresseyi-annotations.gff3
counted the features in the gff (there are no non-genic regions annotated e.g. reapeats):

Feature	num. of features
CDS	30022
exon	30022
gene	23686
mRNA	23686

Checked for features overlapping with CpGs methylated >=50%

Sample	mCpG(>= 50%)	mCpG overlapping with genes/mRNA	mCpG overlapping with exon/CDS
F1	2335	114	106
F2	1864	103	95
F1.merged	1342	60	58
F2.merged	1102	45	42

Moving forward

Overall most mCpGs are not located in genic regions so it’s hard to say how to target this sparse methylation aside from MBD
Would be great to get repeat regions and other features (UTR, etc)
Options for 1 sequencing run:
- resequence 2 individuals to attempt to acheive 100% genome coverage
  - this would give at least 500M reads more data per individual
  - cost: $4,940 (1 Novaseq run)
- WGBS 4 individuals aiming for 400M reads each to attempt acheive 5x coverage of >95% of genome
  - cost: ~$5150 = $4,940 (1 Novaseq run) + library prep (~ $50/sample) + time
  - logic: 226M reads gave 65% genome coverage @5x read depth, 168M gave 45% genome coverage @5x read depth
    - assuming a linear relationship between read depth and genome coverage (2.9 * 100) + 37.5 = ~340M ; see chart below
- WGBS 2-3 individuals aiming for 500M reads each to attempt to achieve 10x coverage of > 95% of genome
  - cost: ~$5100 = $4,940 (1 Novaseq run) + library prep (~ $50/sample) + time
  - logic: 226M reads gave 46% genome coverage @ 10x read depth, 168M gave 31% genome coverage @10x read depth
    - assuming a linear relationship between read depth and genome coverage (3.77 * 100) + 51.8 = ~430M ; see chart below
- MBD-BS on 10 individuals from each population:
  - cost: ~$6600 = $4940 (1 Novaseq run) + library prep (~ $50/sample) + methylminer prep (~ $32/sample) + time
  - could follow optimizations used here: https://www.tandfonline.com/doi/full/10.1080/15592294.2017.1335849

from shellytrigg https://ift.tt/2KKAxbM
via IFTTT

genefish 23:38 on April 30

Grace’s Notebook: Annotate 2019 infection status DEG list; short-term paper plan

Today I read a lot about crustacean immune systems and macropinocytosis. Notes in my crab paper. Additionally, SR and I met and chatted about short-term plan. I annotated the 2019 infection status DEG list. Details in post

Annotating DEG list from 2019 infection status comparison

Rmd: annotate-2019infection_status-DEGs.Rmd

Result: 2019-infection_annot-DEGlist.tab

Need for the file was pointed out by Steven when discussing results section of the crab paper. Added link to the file in the paper.

New short-term paper and crab project plan:

Goal is to focus on the paper and smooth it out. Get Intro and other sections in a more paragraph form.

Talk a little about macropinocytosis now, but don’t get hung up on it. Can go back to it later after the big picture things are written and described.

After I get things more smooth in the paper, I’ll look for some similar papers in crustacea or crab and compare the transcriptome assemblies by number of genes ID-ed, average size and N50s, etc.

FOCUSING ON BIG PICUTRE PAPER PROGRESS CURRENTLY.

from Grace’s Lab Notebook https://ift.tt/2YnAzOW
via IFTTT

genefish 19:24 on April 30

CHiP-Seq and ATAC-Seq

https://bioinformatics-core-shared-training.github.io/cruk-autumn-school-2017/ChIP/Materials/Lectures/Lecture4_Introduction%20to%20ChIP-seq%20and%20ATAC-seq_SS.pdf

genefish 16:26 on April 30

Sam’s Notebook: Transcript Abundance – C.bairdi Alignment-free with Salmon Using 2020-GW Data on Mox

Clarified with Steven an approach for tackling multi-condition comparisons (see this GitHub Issue). As such, I need to have individual transcript abundances for each sample from the 2020 Genewiz RNAseq data before I can proceed. So, I ran salmon (v1.2.1) to perform an alignment-free set of transcript abundances. It’s ridiculously fast, btw…

This was run on Mox and used the C.bairdi-specific reads that were extracted using MEGAN6 on 202020330.

SBATCH script (GitHub):

20200429_cbai_salmon_2020GW_transcript_abundances.sh

#!/bin/bash ## Job Name #SBATCH --job-name=cbai_DEG_basic ## Allocation Definition #SBATCH --account=coenv #SBATCH --partition=coenv ## Resources ## Nodes #SBATCH --nodes=1 ## Walltime (days-hours:minutes:seconds format) #SBATCH --time=10-00:00:00 ## Memory per node #SBATCH --mem=120G ##turn on e-mail notification #SBATCH --mail-type=ALL #SBATCH --mail-user=samwhite@uw.edu ## Specify the working directory for this job #SBATCH --chdir=/gscratch/scrubbed/samwhite/outputs/20200429_cbai_salmon_2020GW_transcript_abundances # Script to generate set of transcript abundances for all C.bairdi Genewiz 2020 data. # # C.bairdi-specific reads were extracted with MEGAN6: # https://robertslab.github.io/sams-notebook/2020/03/30/RNAseq-Reads-Extractions-C.bairdi-Taxonomic-Reads-Extractions-with-MEGAN6-on-swoose.html # # Transcriptome was produced here: https://robertslab.github.io/sams-notebook/2020/03/30/Transcriptome-Assembly-C.bairdi-with-MEGAN6-Taxonomy-specific-Reads-with-Trinity-on-Mox.html # Transcriptome is the same as: cbai_transcriptome_v1.5.fasta # # Salmon index generated during a previous gene expression analysis: # https://robertslab.github.io/sams-notebook/2020/04/22/Gene-Expression-C.bairdi-Pairwise-DEG-Comparisons-with-2019-RNAseq-using-Trinity-Salmon-EdgeR-on-Mox.html ################################################################################################################### # BEGIN USER SETTINGS # Programs array declare -A programs_array programs_array=([salmon]="/gscratch/srlab/programs/salmon-1.2.1_linux_x86_64/bin/salmon") ## Designate input files and locations fastq_dir="/gscratch/srlab/sam/data/C_bairdi/RNAseq/" salmon_index="/gscratch/srlab/sam/data/C_bairdi/transcriptomes/20200408.C_bairdi.megan.Trinity.fasta.salmon.idx" # Set number of CPU threads # Salmon default is 56 threads - so not needed # threads=28 # END USER SETTINGS #################################################################################################################### # Exit script if any command fails set -e # Load Python Mox module for Python module availability module load intel-python3_2017 # Document programs in PATH (primarily for program version ID) { date echo "" echo "System PATH for $SLURM_JOB_ID" echo "" printf "%0.s-" {1..10} echo "${PATH}" | tr : \\n } >> system_path.log # Caputure working directory #wd="$(pwd)" # Capture program options ## NOTE: This particular instance is specific to salmon! for program in "${!programs_array[@]}" do { echo "Program options for ${programs_array[$program]}: " echo "" ${programs_array[$program]} quant --help echo "" ${programs_array[$program]} quant --help-reads echo "" echo "

genefish 01:38 on April 30

Grace’s Notebook: Made new crab transcriptome stress response gene list; notes on crab immunity and macropinocytosis

Today I made a new crab stress response table with more information. I also started an outline for the intro based on the meeting I had with Steven yesterday. I’ve been reading/finding resources about crab/crustacean immunity and macropinocytosis. Will add more to intro and discussion as I learn more.

New stress response gene list

Rmd: 042720-get_crab-stress_response-genelist.Rmd

New list: crab_stress-response-genes.tab

The new list now has BLAST output from the crab transcriptome and the Uniprot-SP-GO information join-ed with the “stress response” list.

Crab paper

Working on:

Intro (reading about crab immunity to add to intro)
Results (will describe transcriptome more)
Discussion (will add more about macropinocytosis and what it means)

from Grace’s Lab Notebook https://ift.tt/2y46RDQ
via IFTTT

genefish 11:57 on April 28

Sam’s Notebook: Go To Goslim C.bairdi Enriched Go Terms From 20200422 Degs

pu— layout: post title: GO to GOslim – C.bairdi Enriched GO Terms from 20200422 DEGs date: ‘2020-04-27 12:04’ tags:

Tanner crab
GO
GOslim
gene ontology
R
Chionoecetes bairdi categories:
Miscellaneous

After running pairwise comparisons and identify differentially expressed genes (DEGs) on 20200422 and finding enriched gene ontology terms, I decided to map the GO terms to Biological Process GOslims. Additionally, I decided to try another level of comparison (I’m not sure how valid it is), whereby I will count the number of GO terms assigned to each GOslim and then calculate the percentage of GOterms that get assigned to each of the GOslim categories. The idea being that it might help identify Biological Processes that are “favored” in a given set of DEGs. I decided to set up “fancy” pyramid plots to view a given set of GO-GOslims for each DEG comparison.

All work was done in R. The initial counting/percentage calculations were done with the following R script (note: all of the followin R code are part of an R Project – link is in the RESULTS section of notebook). The R script uses a “flattened” set of the enriched GO terms identified by Trinity/GOseq, where flattened means one GO term per row. So, a gene may be represented multiple times, in multiple rows if there were multiple GO terms assigned to it by Trinity/GOseq.

The script then relies on the GSEABase package (Bioconductor) and the GO Consortium’s “Generic GO subset”:

goslim_generic.obo

After that, I plotted the outputs.

GO to GOslim R script:

GO_to_GOslim.R

library(GSEABase) library(tidyverse) ######################################################################### # Scipt to map C.bairdi DEG enriched GO terms to GOslims # and identify the GO terms contributing to each GOslim ######################################################################### ### Download files and specify destination directory download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/D9-D12/edgeR.169728.dir/salmon.gene.counts.matrix.D12_vs_D9.edgeR.DE_results.P0.05_C1.D12-UP.subset.GOseq.enriched.flattened", destfile = "./data/D9-D12/salmon.gene.counts.matrix.D12_vs_D9.edgeR.DE_results.P0.05_C1.D12-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/D9-D12/edgeR.169728.dir/salmon.gene.counts.matrix.D12_vs_D9.edgeR.DE_results.P0.05_C1.D9-UP.subset.GOseq.enriched.flattened", destfile = "./data/D9-D12/salmon.gene.counts.matrix.D12_vs_D9.edgeR.DE_results.P0.05_C1.D9-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/D9-D26/edgeR.200352.dir/salmon.gene.counts.matrix.D26_vs_D9.edgeR.DE_results.P0.05_C1.D9-UP.subset.GOseq.enriched.flattened", destfile = "./data/D9-D26/salmon.gene.counts.matrix.D26_vs_D9.edgeR.DE_results.P0.05_C1.D9-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/D9-D26/edgeR.200352.dir/salmon.gene.counts.matrix.D26_vs_D9.edgeR.DE_results.P0.05_C1.D26-UP.subset.GOseq.enriched.flattened", destfile = "./data/D9-D26/salmon.gene.counts.matrix.D26_vs_D9.edgeR.DE_results.P0.05_C1.D26-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/D12-D26/edgeR.230922.dir/salmon.gene.counts.matrix.D12_vs_D26.edgeR.DE_results.P0.05_C1.D26-UP.subset.GOseq.enriched.flattened", destfile = "./data/D12-D26/salmon.gene.counts.matrix.D12_vs_D26.edgeR.DE_results.P0.05_C1.D26-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/ambient-cold/edgeR.267393.dir/salmon.gene.counts.matrix.ambient_vs_cold.edgeR.DE_results.P0.05_C1.ambient-UP.subset.GOseq.enriched.flattened", destfile = "./data/ambient-cold/salmon.gene.counts.matrix.ambient_vs_cold.edgeR.DE_results.P0.05_C1.ambient-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/ambient-cold/edgeR.267393.dir/salmon.gene.counts.matrix.ambient_vs_cold.edgeR.DE_results.P0.05_C1.cold-UP.subset.GOseq.enriched.flattened", destfile = "./data/ambient-cold/salmon.gene.counts.matrix.ambient_vs_cold.edgeR.DE_results.P0.05_C1.cold-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/ambient-warm/edgeR.297991.dir/salmon.gene.counts.matrix.ambient_vs_warm.edgeR.DE_results.P0.05_C1.warm-UP.subset.GOseq.enriched.flattened", destfile = "./data/ambient-warm/salmon.gene.counts.matrix.ambient_vs_warm.edgeR.DE_results.P0.05_C1.warm-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/ambient-warm/edgeR.297991.dir/salmon.gene.counts.matrix.ambient_vs_warm.edgeR.DE_results.P0.05_C1.ambient-UP.subset.GOseq.enriched.flattened", destfile = "./data/ambient-warm/salmon.gene.counts.matrix.ambient_vs_warm.edgeR.DE_results.P0.05_C1.ambient-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/cold-warm/edgeR.328585.dir/salmon.gene.counts.matrix.cold_vs_warm.edgeR.DE_results.P0.05_C1.warm-UP.subset.GOseq.enriched.flattened", destfile = "./data/cold-warm/salmon.gene.counts.matrix.cold_vs_warm.edgeR.DE_results.P0.05_C1.warm-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/cold-warm/edgeR.328585.dir/salmon.gene.counts.matrix.cold_vs_warm.edgeR.DE_results.P0.05_C1.cold-UP.subset.GOseq.enriched.flattened", destfile = "./data/cold-warm/salmon.gene.counts.matrix.cold_vs_warm.edgeR.DE_results.P0.05_C1.cold-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/infected-uninfected/edgeR.132470.dir/salmon.gene.counts.matrix.infected_vs_uninfected.edgeR.DE_results.P0.05_C1.infected-UP.subset.GOseq.enriched.flattened", destfile = "./data/infected-uninfected/salmon.gene.counts.matrix.infected_vs_uninfected.edgeR.DE_results.P0.05_C1.infected-UP.subset.GOseq.enriched.flattened") download.file(url = "https://gannet.fish.washington.edu/Atumefaciens/20200422_cbai_DEG_basic_comparisons/infected-uninfected/edgeR.132470.dir/salmon.gene.counts.matrix.infected_vs_uninfected.edgeR.DE_results.P0.05_C1.uninfected-UP.subset.GOseq.enriched.flattened", destfile = "./data/infected-uninfected/salmon.gene.counts.matrix.infected_vs_uninfected.edgeR.DE_results.P0.05_C1.uninfected-UP.subset.GOseq.enriched.flattened") ### Set false discovery rate (FDR) filter, if desired fdr <- as.character("1.0") ### Create list of files goseq_files <- list.files(path = "./data", pattern = "\\.GOseq.[de]", recursive = TRUE, full.names = TRUE) ### Set output filename suffix output_suffix=("GOslims.csv") ### Strip path from goseq files goseq_filename <- basename(goseq_files) ### Vector of GOslim ontologies (e.g. Biological Process = BP, Molecular Function = MF, Cellular Component = CC) ontologies <- c("BP", "CC", "MF") for (slim_ontology in ontologies) { ### Set GOOFFSPRING database, based on ontology group set above go_offspring <- paste("GO", slim_ontology, "OFFSPRING", sep = "") for (item in goseq_files) { ## Get max number of fields # Needed to handle reading in file with different number of columns in each row # as there may be multiple max_fields <- max(count.fields(item, sep = "\t"), na.rm = TRUE) ## Read in tab-delimited GOseq file # Use "max_fields" to populate all columns with a sequentially numbered header go_seqs <- read.delim(item, col.names = paste0("V",seq_len(max_fields))) ## Filter enriched GOterms with false discovery rate goseqs_fdr <- filter(go_seqs, V8 <= as.numeric(fdr)) ## Grab just the individual GO terms from the "category" column) goterms <- as.character(goseqs_fdr$V1) ### Use GSEA to map GO terms to GOslims ## Store goterms as GSEA object myCollection <- GOCollection(goterms) ## Use generic GOslim file to create a GOslim collection # I downloaded goslim_generic.obo from http://geneontology.org/docs/go-subset-guide/ # then i moved it to the R library for GSEABase in the extdata folder # in addition to using the command here - I think they're both required. slim <- getOBOCollection("./data/goslim_generic.obo") ## Map GO terms to GOslims and select Biological Processes group slimsdf <- goSlim(myCollection, slim, slim_ontology) ## Need to know the 'offspring' of each term in the ontology, and this is given by the data in: # GO.db::getFromNamespace(go_offspring, "GO.db") ## Create function to parse out GO terms assigned to each GOslim ## Courtesy Bioconductor Support: https://support.bioconductor.org/p/128407/ mappedIds <- function(df, collection, OFFSPRING) { map <- as.list(OFFSPRING[rownames(df)]) mapped <- lapply(map, intersect, ids(collection)) df[["go_terms"]] <- vapply(unname(mapped), paste, collapse = ";", character(1L)) df } ## Run the function slimsdf <- mappedIds(slimsdf, myCollection, getFromNamespace(go_offspring, "GO.db")) ## Provide column name for first column slimsdf <- cbind(GOslim = rownames(slimsdf), slimsdf) rownames(slimsdf) <- NULL ### Prep output file naming structure ## Split filenames on periods ## Creates a list split_filename <- strsplit(item, ".", fixed = TRUE) ## Split filename on directories ## Creates a list split_dirs <- strsplit(item, "/", fixed = TRUE) ## Slice split_filename list from position 9 to the last position of the list ## Paste these together using a period goseq_filename_split <-paste(split_filename[[1]][9:lengths(split_filename)], collapse = ".") ## Slice split_dirs list at position ## Paste elements together to form output filename fdr_file_out <- paste("FDR", fdr, sep = "_") outfilename <- paste(goseq_filename_split, fdr_file_out, slim_ontology ,output_suffix, collapse = ".", sep = ".") ## Set output file destination and name ## Adds proper subdirectory from split_dirs list outfile_dest <- file.path("./analyses", split_dirs[[1]][3], outfilename) ## Write output file write.csv(slimsdf, file = outfile_dest, quote = FALSE, row.names = FALSE) } }

Since I got “lazy” and didn’t want to try to figure out how to properly loop through all of the output files from the above script, I just made individual scripts to plot each set of comparison GO terms percentages assigned to GOslims. Here’s an example:

infected-uninfected_GOslim_pyramid_plotting.R

# Script to generate a "pyramid" plot # comparing the percentages of enriched GO terms assinged # to each category of Biological Process GOslims library(dplyr) library(ggplot2) ##################################################### # Set the following variables for the appropriate comparisons/files ##################################################### # Comparison comparison <- "infected-uninfected" # Treatments treatment_01 <- "infected" treatment_02 <- "uninfected" # Read in first comparsion files df1 <- read.csv("analyses/infected-uninfected/P0.05_C1.infected-UP.subset.GOseq.enriched.flattened.FDR_1.0.BP.GOslims.csv") # Read in second comparison file df2 <- read.csv("analyses/infected-uninfected/P0.05_C1.uninfected-UP.subset.GOseq.enriched.flattened.FDR_1.0.BP.GOslims.csv") ###################################################### # CHANGES BELOW HERE ARE PROBABLY NOT NECESSARY ###################################################### # GOslim categories ontologies <- c("BP", "CC", "MF") # Remove generic "biological_process" category df1 <- df1[df1$GOslim != "GO:0008150",] df2 <- df2[df2$GOslim != "GO:0008150",] # Remove generic "cellular_component" category df1 <- df1[df1$GOslim != "GO:0005575",] df2 <- df2[df2$GOslim != "GO:0005575",] # Remove generic "molecular_function" category df1 <- df1[df1$GOslim != "GO:0003674",] df2 <- df2[df2$GOslim != "GO:0003674",] # Select columns df1 <- df1 %>% select(Term, Percent) df2 <- df2 %>% select(Term, Percent) # Create treatment column and assign term to all rows df1$treatment <- treatment_01 df2$treatment <- treatment_02 # Concatenate dataframes df3 <- rbind(df1, df2) # Filename for plot pyramid <- paste(comparison, "GOslims", "BP", "png", sep = ".") pyramid_path <- paste(comparison, pyramid, sep = "/") pyramid_dest <- file.path("./analyses", pyramid_path) # "Open" PNG file for saving subsequent plot png(pyramid_dest, width = 600, height = 1200, units = "px", pointsize = 12) # Create "pyramid" plot ggplot(df3, aes(x = Term, fill = treatment, y = ifelse(test = treatment == treatment_01, yes = -Percent, no = Percent))) + geom_bar(stat = "identity") + scale_y_continuous(labels = abs, limits = max(df3$Percent) * c(-1,1)) + labs(title = "Percentages of GO terms assigned to BP GOslims", x = "GOslim", y = "Percent GO terms in GOslim") + scale_x_discrete(expand = c(-1,0)) + coord_flip() # Close PNG file dev.off()

genefish 05:59 on April 28

Yaamini’s Notebook: MethCompare Part 8

Characterizing CpGs in 5x data

After confirming last week that mapping to a pan-genome wasn’t changing our output, I went ahead with the CpG characterization pipeline. We decided to move ahead with 5x data with the option to return to 10x data later. I decided to create a new Jupyter notebook for running 5x data through the pipeline.

Repo fiasco

Turns out we haven’t mastered how to have five people working on a repo simultaneously without having pull-and-push issues. When Sam reverted to an old repo version, I lost the checksum code I worked on earlier. The first thing I did was add the code back in to my 10x Jupyter notebook and to the new notebook I created for 5x data. This isn’t the first time a repo commit from someone else lost some of my code, so hopefully we figure out!

M. capitata

Using 5x files provided more CpG data to characterize, but the percentages of methylated CpGs total and in various genome features were mostly consistent between 5x and 10x data.

Table 1. Methylated, sparsely methylated, and unmethylated CpGs in M. capitata

Sample	Method	CpG with Data	Methylated CpG	Sparsely Methylated CpG	Unmethylated CpG
10	WGBS	4571288	450582 (9.9%)	547868 (12.0%)	3572838 (78.2%)
11	WGBS	4661716	528902 (11.3%)	517805 (11.1%)	3615009 (77.5%)
12	WGBS	8791700	1059904 (12.1%)	1000337 (11.4%)	6731459 (76.6%)
13	RRBS	3173254	257741 (8.1%)	152042 (4.8%)	2763471 (87.1%)
14	RRBS	2648697	184742 (7.0%)	135052 (5.1%)	2328903 (87.9%)
15	RRBS	3176517	231347 (7.3%)	179454 (5.6%)	2765716 (87.1%)
16	MBD-BSSeq	583599	106695 (18.3%)	74839 (12.8%)	402065 (68.9%)
17	MBD-BSSeq	242390	45506 (18.8%)	28850 (11.9%)	168034 (69.3%)
18	MBD-BSSeq	153392	29468 (19.2%)	16793 (10.9%)	107131 (69.8%)

Table 2. Number and percent of CpGs that overlap with genes in M. capitata.

Sample	Method	CpGs with Data	Methylated CpGs	Sparsely Methylated CpGs	Unmethylated CpGs
10	WGBS	1683069	230343 (13.7%)	196827 (11.7%)	1255899 (74.6%)
11	WGBS	1740149	267181 (15.3%)	188348 (10.8%)	1284620 (73.8%)
12	WGBS	3204885	533230 (16.6%)	348967 (10.9%)	2322688 (72.5%)
13	RRBS	1062577	123271 (11.6%)	53235 (5.0%)	886071 (83.4%)
14	RRBS	895852	89865 (10.0%)	48499 (5.4%)	757488 (84.6%)
15	RRBS	1064769	113217 (10.6%)	64106 (6.0%)	887446 (83.3%)
16	MBD-BSSeq	208503	48436 (23.2%)	25316 (12.1%)	134751 (64.6%)
17	MBD-BSSeq	83519	20458 (24.5%)	9441 (11.3%)	53620 (64.2%)
18	MBD-BSSeq	50302	12960 (25.8%)	4843 (9.6%)	32499 (64.6%)

Table 3. Number and percent of CpGs that overlap with CDS in M. capitata.

Sample	Method	CpGs with Data	Methylated CpGs	Sparsely Methylated CpGs	Unmethylated CpGs
10	WGBS	476579	54412 (11.4%)	60266 (12.6%)	361901 (75.9%)
11	WGBS	496887	64070 (12.9%)	58258 (11.7%)	374559 (75.4%)
12	WGBS	791965	113396 (14.3%)	89455 (11.3%)	589114 (74.4%)
13	RRBS	200808	18158 (9.0%)	11383 (5.7%)	171267 (85.3%)
14	RRBS	172313	12833 (7.4%)	10486 (6.1%)	148994 (86.5%)
15	RRBS	203147	15130 (7.4%)	14015 (6.9%)	174002 (85.7%)
16	MBD-BSSeq	72490	13535 (18.7%)	9666 (13.3%)	49289 (68.0%)
17	MBD-BSSeq	29305	5649 (19.3%)	3690 (12.6%)	19966 (68.1%)
18	MBD-BSSeq	17619	3992 (22.7%)	1907 (10.8%)	11720 (66.5%)

Table 4. Number and percent of CpGs that overlap with introns in M. capitata.

Sample	Method	CpGs with Data	Methylated CpGs	Sparsely Methylated CpGs	Unmethylated CpGs
10	WGBS	1207609	176102 (14.6%)	136671 (11.3%)	894836 (74.1%)
11	WGBS	1244403	203311 (16.3%)	130195 (10.5%)	910897 (73.2%)
12	WGBS	2415036	420249 (17.4%)	259702 (10.8%)	1735085 (71.8%)
13	RRBS	862191	105145 (12.2%)	41866 (4.9%)	715180 (82.9%)
14	RRBS	723933	77053 (10.6%)	38032 (5.3%)	608848 (84.1%)
15	RRBS	862091	98111 (11.4%)	50116 (5.8%)	713864 (82.8%)
16	MBD-BSSeq	136153	34929 (25.7%)	15666 (11.5%)	85558 (62.8%)
17	MBD-BSSeq	54275	14821 (27.3%)	5755 (10.6%)	33699 (62.1%)
18	MBD-BSSeq	32717	8973 (27.4%)	2942 (9.0%)	20802 (63.6%)

Table 5. Number and percent of CpGs that overlap with intergenic regions in M. capitata.

Sample	Method	CpGs with Data	Methylated CpGs	Sparsely Methylated CpGs	Unmethylated CpGs
10	WGBS	2888219	220239 (7.6%)	351041 (12.2%)	2316939 (80.2%)
11	WGBS	2921567	261721 (9.0%)	329457 (11.3%)	2330389 (79.8%)
12	WGBS	5586815	526674 (9.4%)	651370 (11.7%)	4408771 (78.9%)
13	RRBS	2110677	134470 (6.4%)	98807 (4.7%)	1877400 (88.9%)
14	RRBS	1752845	94877 (5.4%)	86553 (4.9%)	1571415 (89.6%)
15	RRBS	2111748	118130 (5.6%)	115348 (5.5%)	1878270 (88.9%)
16	MBD-BSSeq	375096	58259 (15.5%)	49523 (13.2%)	267314 (71.3%)
17	MBD-BSSeq	158871	25048 (15.8%)	19409 (12.2%)	114414 (72.0%)
18	MBD-BSSeq	103090	16508 (16.0%)	11950 (11.6%)	74632 (72.4%)

P. acuta

Table 6. Methylated, sparsely methylated, and unmethylated CpGs in P. acuta

Sample	Method	CpG with Data	Methylated CpG	Sparsely Methylated CpG	Unmethylated CpG
1	WGBS	5546051	110364 (2.0%)	367019 (6.6%)	5068668 (91.4%)
2	WGBS	6358722	126440 (2.0%)	345887 (5.4%)	5886395 (92.6%)
3	WGBS	5866786	124819 (2.1%)	385346 (6.6%)	5356621 (91.3%)
4	RRBS	1835561	31047 (1.7%)	137700 (7.5%)	1666814 (90.8%)
5	RRBS	1451229	30345 (2.1%)	64837 (4.5%)	1356047 (93.4%)
6	RRBS	1517358	26617 (1.8%)	89246 (5.9%)	1401495 (92.4%)
7	MBD-BSSeq	2640625	258222 (9.8%)	296059 (11.2%)	2086344 (79.0%)
8	MBD-BSSeq	539008	213342 (39.6%)	80086 (14.9%)	245580 (45.6%)
9	MBD-BSSeq	2732607	255370 (9.3%)	337855 (12.4%)	2139382 (78.3%)

Table 7. Number and percent of CpGs that overlap with genes in P. acuta.

Sample	Method	CpG with Data	Methylated CpG	Sparsely Methylated CpG	Unmethylated CpG
1	WGBS	2466992	73959 (3.0%)	157337 (6.4%)	2235696 (90.6%)
2	WGBS	2761956	85861 (3.1%)	144292 (5.2%)	2531803 (91.7%)
3	WGBS	2588278	82377 (3.2%)	161791 (6.3%)	2344110 (90.6%)
4	RRBS	776123	13588 (1.8%)	56290 (7.3%)	706245 (91.0%)
5	RRBS	607225	12789 (2.1%)	25954 (4.3%)	568482 (93.6%)
6	RRBS	639570	11396 (1.8%)	35633 (5.6%)	592541 (92.6%)
7	MBD-BSSeq	1253805	118291 (9.4%)	120485 (9.6%)	1015029 (81.0%)
8	MBD-BSSeq	220096	86074 (39.1%)	27902 (12.7%)	106120 (48.2%)
9	MBD-BSSeq	1281644	125536 (9.8%)	139034 (10.8%)	1017074 (79.4%)

Table 8. Number and percent of CpGs that overlap with CDS in P. acuta.

Sample	Method	CpG with Data	Methylated CpG	Sparsely Methylated CpG	Unmethylated CpG
1	WGBS	1494340	59188 (4.0%)	89863 (6.0%)	1345289 (90.0%)
2	WGBS	1620632	66365 (4.1%)	76868 (4.7%)	1477399 (91.2%)
3	WGBS	1552715	65245 (4.2%)	89654 (5.8%)	1397816 (90.0%)
4	RRBS	500779	9644 (1.9%)	36616 (7.3%)	454519 (90.8%)
5	RRBS	395869	8808 (2.2%)	17345 (4.4%)	369716 (93.4%)
6	RRBS	415529	7954 (1.9%)	23679 (5.7%)	383896 (92.4%)
7	MBD-BSSeq	931250	92559 (9.9%)	87369 (9.4%)	751322 (80.7%)
8	MBD-BSSeq	184606	70696 (38.3%)	21973 (11.9%)	91937 (49.8%)
9	MBD-BSSeq	924825	95614 (10.3%)	98699 (10.7%)	730512 (79.0%)

Table 9. Number and percent of CpGs that overlap with introns in P. acuta.

Sample	Method	CpG with Data	Methylated CpG	Sparsely Methylated CpG	Unmethylated CpG
1	WGBS	1840138	41787 (2.3%)	122271 (6.6%)	1676080 (91.1%)
2	WGBS	2113295	51567 (2.4%)	117738 (5.6%)	1943990 (92.0%)
3	WGBS	1946150	47352 (2.4%)	128122 (6.6%)	1770676 (91.0%)
4	RRBS	539524	8446 (1.6%)	38589 (7.2%)	492489 (91.3%)
5	RRBS	414147	8375 (2.0%)	17239 (4.2%)	388533 (93.8%)
6	RRBS	439844	7162 (1.6%)	23825 (5.4%)	408857 (93.0%)
7	MBD-BSSeq	746710	64572 (8.6%)	71046 (9.5%)	611092 (81.8%)
8	MBD-BSSeq	101087	42068 (41.6%)	13365 (13.2%)	45654 (45.2%)
9	MBD-BSSeq	794295	72530 (9.1%)	85117 (10.7%)	636648 (80.2%)

Table 10. Number and percent of CpGs that overlap with intergenic regions in P. acuta.

Sample	Method	CpG with Data	Methylated CpG	Sparsely Methylated CpG	Unmethylated CpG
1	WGBS	3080835	3080835 (1.2%)	209781 (6.8%)	2834593 (92.0%)
2	WGBS	3598848	40642 (1.1%)	201712 (5.6%)	3356494 (93.3%)
3	WGBS	3280357	42507 (1.3%)	223666 (6.8%)	3014184 (91.9%)
4	RRBS	1060062	17473 (1.6%)	81473 (7.7%)	961116 (90.7%)
5	RRBS	844509	17581 (2.1%)	38900 (4.6%)	788028 (93.3%)
6	RRBS	878266	15232 (1.7%)	53637 (6.1%)	809397 (92.2%)
7	MBD-BSSeq	1387623	140057 (10.1%)	175668 (12.7%)	1071898 (77.2%)
8	MBD-BSSeq	319125	127376 (39.9%)	52215 (16.4%)	139534 (43.7%)
9	MBD-BSSeq	1451853	129949 (9.0%)	198940 (13.7%)	1122964 (77.3%)

Going forward

Create figures for CpG characterization in various genome features
Update code for methylation frequency distribution figure
Generate repeat tracks for each species
Find a way to generate tables programmatically
Figure out how to meaningfully concatenate data for each method
Figure out methylation island analysis

Please enable JavaScript to view the comments powered by Disqus.

from the responsible grad student https://ift.tt/35ak8GU
via IFTTT

genefish 18:57 on April 27

Sam’s Notebook: Gene Expression – C.bairdi Pairwise DEG Comparisons with 2019 RNAseq using Trinity-Salmon-EdgeR on Mox

Per a Slack request, Steven asked me to take the Genewize RNAseq data (received 2020318) through edgeR. Ran the analysis using the Trinity differential expression pipeline:

Here’re the core input files used for this analysis:

Transcriptome: cbai_transcriptome_v1.5.fasta
MEGAN6 Arthropoda taxonomic reads: 20200413_C_bairdi_megan_reads/

The analyses will perform the following pairwise comparisons:

infected-uninfected
D9-D12
D9-D26
D12-D26
ambient-cold
ambient-warm
cold-warm

It will identify differentially expressed genes with >=2-fold log change in expression and a false discovery rate of <=0.05. Additionally, it will perform gene ontology (GO) enrichment analysis using GOseq.

As a brief aside, I’m pretty stoked about the SBATCH script below! It automates FastQ file selection for each comparison, creates appropriately named subdirectories and creates proper Trinity samples list file needed.

After running the DEG analysis, I “flattened” the enriched GO terms files for later use in R to map these GO terms to GOslims. That was run separately and the script is after the SBATCH script.

SBATCH script (GitHub):

20200422_cbai_DEG_basic_comparisons.sh

#!/bin/bash ## Job Name #SBATCH --job-name=cbai_DEG_basic ## Allocation Definition #SBATCH --account=coenv #SBATCH --partition=coenv ## Resources ## Nodes #SBATCH --nodes=1 ## Walltime (days-hours:minutes:seconds format) #SBATCH --time=10-00:00:00 ## Memory per node #SBATCH --mem=120G ##turn on e-mail notification #SBATCH --mail-type=ALL #SBATCH --mail-user=samwhite@uw.edu ## Specify the working directory for this job #SBATCH --chdir=/gscratch/scrubbed/samwhite/outputs/20200422_cbai_DEG_basic_comparisons # This is a script to identify differentially expressed genes (DEGs) in C.bairdi # using pairwise comparisions of from just the "2020-GW" (i.e. just Genewiz) RNAseq data # which has been taxonomically selected for all Arthropoda reads. See Sam's notebook from 20200419 # https://robertslab.github.io/sams-notebook/ # Script will run Trinity's builtin differential gene expression analysis using: # - Salmon alignment-free transcript abundance estimation # - edgeR # Cutoffs of 2-fold difference in expression and FDR of <=0.05. ################################################################################### # These variables need to be set by user fastq_dir="/gscratch/srlab/sam/data/C_bairdi/RNAseq/" fasta_prefix="20200408.C_bairdi.megan.Trinity" transcriptome_dir="/gscratch/srlab/sam/data/C_bairdi/transcriptomes" trinotate_feature_map="${transcriptome_dir}/20200409.cbai.trinotate.annotation_feature_map.txt" go_annotations="${transcriptome_dir}/20200409.cbai.trinotate.go_annotations.txt" # Array of the various comparisons to evaluate # Each condition in each comparison should be separated by a "-" comparisons_array=( infected-uninfected \ D9-D12 \ D9-D26 \ D12-D26 \ ambient-cold \ ambient-warm \ cold-warm ) # Functions # Expects input (i.e. "$1") to be in the following format: # e.g. 20200413.C_bairdi.113.D9.uninfected.cold.megan_R2.fq get_day () { day=$(echo "$1" | awk -F"." '{print $4}'); } get_inf () { inf=$(echo "$1" | awk -F"." '{print $5}'); } get_temp () { temp=$(echo "$1" | awk -F"." '{print $6}'); } ################################################################################### # Exit script if any command fails set -e # Load Python Mox module for Python module availability module load intel-python3_2017 # Document programs in PATH (primarily for program version ID) { date echo "" echo "System PATH for $SLURM_JOB_ID" echo "" printf "%0.s-" {1..10} echo "${PATH}" | tr : \\n } >> system_path.log wd="$(pwd)" threads=28 ## Designate input file locations transcriptome="${transcriptome_dir}/${fasta_prefix}.fasta" fasta_seq_lengths="${transcriptome_dir}/${fasta_prefix}.fasta.seq_lens" gene_map="${transcriptome_dir}/${fasta_prefix}.fasta.gene_trans_map" transcriptome="${transcriptome_dir}/${fasta_prefix}.fasta" # Standard output/error files diff_expr_stdout="diff_expr_stdout.txt" diff_expr_stderr="diff_expr_stderr.txt" matrix_stdout="matrix_stdout.txt" matrix_stderr="matrix_stderr.txt" salmon_stdout="salmon_stdout.txt" salmon_stderr="salmon_stderr.txt" tpm_length_stdout="tpm_length_stdout.txt" tpm_length_stderr="tpm_length_stderr.txt" trinity_DE_stdout="trinity_DE_stdout.txt" trinity_DE_stderr="trinity_DE_stderr.txt" edgeR_dir="" #programs trinity_home=/gscratch/srlab/programs/trinityrnaseq-v2.9.0 trinity_annotate_matrix="${trinity_home}/Analysis/DifferentialExpression/rename_matrix_feature_identifiers.pl" trinity_abundance=${trinity_home}/util/align_and_estimate_abundance.pl trinity_matrix=${trinity_home}/util/abundance_estimates_to_matrix.pl trinity_DE=${trinity_home}/Analysis/DifferentialExpression/run_DE_analysis.pl diff_expr=${trinity_home}/Analysis/DifferentialExpression/analyze_diff_expr.pl trinity_tpm_length=${trinity_home}/util/misc/TPM_weighted_gene_length.py # Loop through each comparison # Will create comparison-specific direcctories and copy # appropriate FastQ files for each comparison. # After file transfer, will create necessary sample list file for use # by Trinity for running differential gene expression analysis and GO enrichment. for comparison in "${!comparisons_array[@]}" do # Assign variables cond1_count=0 cond2_count=0 comparison=${comparisons_array[${comparison}]} comparison_dir=${wd}/${comparison}/ salmon_gene_matrix=${comparison_dir}/salmon.gene.TMM.EXPR.matrix salmon_iso_matrix=${comparison_dir}/salmon.isoform.TMM.EXPR.matrix samples=${comparison_dir}${comparison}.samples.txt # Extract each comparison from comparisons array # Conditions must be separated by a "-" cond1=$(echo "${comparison}" | awk -F"-" '{print $1}') cond2=$(echo "${comparison}" | awk -F"-" '{print $2}') mkdir --parents "${comparison}" cd "${comparison}" || exit # Series of if statements to identify which FastQ files to rsync to working directory if [[ "${comparison}" == "infected-uninfected" ]]; then rsync --archive --verbose ${fastq_dir}*.fq . fi if [[ "${comparison}" == "D9-D12" ]]; then for fastq in "${fastq_dir}"*.fq do get_day "${fastq}" if [[ "${day}" == "D9" || "${day}" == "D12" ]]; then rsync --archive --verbose "${fastq}" . fi done fi if [[ "${comparison}" == "D9-D26" ]]; then for fastq in "${fastq_dir}"*.fq do get_day "${fastq}" if [[ "${day}" == "D9" || "${day}" == "D26" ]]; then rsync --archive --verbose "${fastq}" . fi done fi if [[ "${comparison}" == "D12-D26" ]]; then for fastq in "${fastq_dir}"*.fq do get_day "${fastq}" if [[ "${day}" == "D12" || "${day}" == "D26" ]]; then rsync --archive --verbose "${fastq}" . fi done fi if [[ "${comparison}" == "ambient-cold" ]]; then #statements for fastq in "${fastq_dir}"*.fq do get_temp "${fastq}" if [[ "${temp}" == "ambient" || "${temp}" == "cold" ]]; then rsync --archive --verbose "${fastq}" . fi done fi if [[ "${comparison}" == "ambient-warm" ]]; then for fastq in "${fastq_dir}"*.fq do get_temp "${fastq}" if [[ "${temp}" == "ambient" || "${temp}" == "warm" ]]; then rsync --archive --verbose "${fastq}" . fi done fi if [[ "${comparison}" == "cold-warm" ]]; then for fastq in "${fastq_dir}"*.fq do get_temp "${fastq}" if [[ "${temp}" == "cold" || "${temp}" == "warm" ]]; then rsync --archive --verbose "${fastq}" . fi done fi # Create reads array # Paired reads files will be sequentially listed in array (e.g. 111_R1 111_R2) reads_array=(*.fq) echo "" echo "Created reads_array" # Loop to create sample list file # Sample file list is tab-delimited like this: # cond_A cond_A_rep1 A_rep1_left.fq A_rep1_right.fq # cond_A cond_A_rep2 A_rep2_left.fq A_rep2_right.fq # cond_B cond_B_rep1 B_rep1_left.fq B_rep1_right.fq # cond_B cond_B_rep2 B_rep2_left.fq B_rep2_right.fq # Increment by 2 to process next pair of FastQ files for (( i=0; i<${#reads_array[@]} ; i+=2 )) do echo "" echo "Evaluating ${reads_array[i]} and ${reads_array[i+1]}" get_day "${reads_array[i]}" get_inf "${reads_array[i]}" get_temp "${reads_array[i]}" echo "" echo "Got day (${day}), infection status (${inf}), and temp (${temp})." echo "" echo "Condition 1 is: ${cond1}" echo "condition 2 is: ${cond2}" # Evaluate specified treatment conditions and format sample file list appropriately. if [[ "${cond1}" == "${day}" || "${cond1}" == "${inf}" || "${cond1}" == "${temp}" ]]; then cond1_count=$((cond1_count+1)) echo "" echo "Condition 1 evaluated." # Create tab-delimited samples file. printf "%s\t%s%02d\t%s\t%s\n" "${cond1}" "${cond1}_" "${cond1_count}" "${comparison_dir}${reads_array[i]}" "${comparison_dir}${reads_array[i+1]}" \ >> "${samples}" elif [[ "${cond2}" == "${day}" || "${cond2}" == "${inf}" || "${cond2}" == "${temp}" ]]; then cond2_count=$((cond2_count+1)) echo "" echo "Condition 2 evaluated." # Create tab-delimited samples file. printf "%s\t%s%02d\t%s\t%s\n" "${cond2}" "${cond2}_" "${cond2_count}" "${comparison_dir}${reads_array[i]}" "${comparison_dir}${reads_array[i+1]}" \ >> "${samples}" fi # Copy sample list file to transcriptome directory cp "${samples}" "${transcriptome_dir}" done echo "Created ${comparison} sample list file." # Create directory/sample list for ${trinity_matrix} command trin_matrix_list=$(awk '{printf "%s%s", $2, "/quant.sf " }' "${samples}") # Determine transcript abundances using Salmon alignment-free # abundance estimate. ${trinity_abundance} \ --output_dir "${comparison_dir}" \ --transcripts ${transcriptome} \ --seqType fq \ --samples_file "${samples}" \ --est_method salmon \ --aln_method bowtie2 \ --gene_trans_map "${gene_map}" \ --prep_reference \ --thread_count "${threads}" \ 1> "${comparison_dir}"${salmon_stdout} \ 2> "${comparison_dir}"${salmon_stderr} # Convert abundance estimates to matrix ${trinity_matrix} \ --est_method salmon \ --gene_trans_map ${gene_map} \ --out_prefix salmon \ --name_sample_by_basedir \ ${trin_matrix_list} \ 1> ${matrix_stdout} \ 2> ${matrix_stderr} # Integrate Trinotate functional annotations "${trinity_annotate_matrix}" \ "${trinotate_feature_map}" \ salmon.gene.counts.matrix \ > salmon.gene.counts.annotated.matrix # Generate weighted gene lengths "${trinity_tpm_length}" \ --gene_trans_map "${gene_map}" \ --trans_lengths "${fasta_seq_lengths}" \ --TPM_matrix "${salmon_iso_matrix}" \ > Trinity.gene_lengths.txt \ 2> ${tpm_length_stderr} # Differential expression analysis # Utilizes edgeR. # Needs to be run in same directory as transcriptome. cd ${transcriptome_dir} || exit ${trinity_DE} \ --matrix "${comparison_dir}salmon.gene.counts.matrix" \ --method edgeR \ --samples_file "${samples}" \ 1> ${trinity_DE_stdout} \ 2> ${trinity_DE_stderr} mv edgeR* "${comparison_dir}" # Run differential expression on edgeR output matrix # Set fold difference to 2-fold (ie. -C 1 = 2^1) # P value <= 0.05 # Has to run from edgeR output directory # Pulls edgeR directory name and removes leading ./ in find output # Using find is required because edgeR names directory using PID # and I don't know how to find that out cd "${comparison_dir}" || exit edgeR_dir=$(find . -type d -name "edgeR*" | sed 's%./%%') cd "${edgeR_dir}" || exit mv "${transcriptome_dir}/${trinity_DE_stdout}" . mv "${transcriptome_dir}/${trinity_DE_stderr}" . ${diff_expr} \ --matrix "${salmon_gene_matrix}" \ --samples "${samples}" \ --examine_GO_enrichment \ --GO_annots "${go_annotations}" \ --include_GOplot \ --gene_lengths "${comparison_dir}Trinity.gene_lengths.txt" \ -C 1 \ -P 0.05 \ 1> ${diff_expr_stdout} \ 2> ${diff_expr_stderr} cd "${wd}" || exit done

Flatten enriched GO terms file (GitHub):

trinity_deg_to_go.sh

#!/bin/bash ############################################################# # Script to "flatten" Trinity edgeR GOseq enrichment format # so each line contains a single gene/transcript ID # and associated GO term ############################################################# # Enable globstar for recursive searching shopt -s globstar # Declare variables output_file="" wd=$(pwd) # Input file ## Expects Trinity edgeR GOseq enrichment format: ## category over_represented_pvalue under_represented_pvalue numDEInCat numInCat term ontology over_represented_FDR go_term gene_ids ## Field 10 (gene_ids) contains comma separated gene_ids that fall in the given GO term in the "category" column for goseq in **/*UP.subset*.enriched do # Capture path to file dir=${goseq%/*} cd "${dir}" || exit tmp_file=$(mktemp) # Count lines in file linecount=$(cat "${goseq}" | wc -l) # If file is not empty if (( "${linecount}" > 1 )) then output_file="${goseq}.flattened" # 1st: Convert comma-delimited gene IDs in column 10 to tab-delimited # Also, set output (OFS) to be tab-delimited # 2nd: Convert spaces to underscores and keep output as tab-delimited # 3rd: Sort on Trinity IDs (column 10) and keep only uniques awk 'BEGIN{FS="\t";OFS="\t"} {gsub(/, /, "\t", $10); print}' "${goseq}" \ | awk 'BEGIN{F="\t";OFS="\t"} NR==1; NR > 1 {gsub(/ /, "_", $0); print}' \ > "${tmp_file}" # Identify the first line number which contains a gene_id begin_goterms=$(grep --line-number "TRINITY" "${tmp_file}" \ | awk '{for (i=1;i<=NF;i++) if($i ~/TRINITY/) print i}' \ | sort --general-numeric-sort --unique | head -n1) # "Unfolds" gene_ids to a single gene_id per row while read -r line do # Capture the length of the longest row max_field=$(echo "$line" | awk -F "\t" '{print NF}') # Retain the first 8 fields (i.e. categories) fixed_fields=$(echo "$line" | cut -f1-8) # Since not all the lines contain the same number of fields (e.g. may not have GO terms), # evaluate the number of fields in each line to determine how to handle current line. # If the value in max_field is less than the field number where the GO terms begin, # then just print the current line (%s) followed by a newline (\n). if (( "$max_field" < "$begin_goterms" )) then printf "%s\n" "$line" else goterms=$(echo "$line" | cut -f"$begin_goterms"-"$max_field") # Assign values in the variable "goterms" to a new indexed array (called "array"), # with tab delimiter (IFS=$'\t') IFS=$'\t' read -r -a array <<<"$goterms" # Iterate through each element of the array. # Print the first n fields (i.e. the fields stored in "fixed_fields") followed by a tab (%s\t). # Print the current element in the array (i.e. the current GO term) followed by a new line (%s\n). for element in "${!array[@]}" do printf "%s\t%s\n" "$fixed_fields" "${array[$element]}" done fi done < "${tmp_file}" > "${output_file}" fi # Cleanup rm "${tmp_file}" cd "${wd}" || exit done

genefish 11:57 on April 27

Sam’s Notebook: FastQC-MultiQC – Laura Spencer’s QuantSeq Data

Laura Spencer received her O.lurida QuantSeq data, so I put it through FastQC/MultiQC and put the pertinent info in the nightingales Google Sheet. I also moved the data to /owl/nightingales/O_lurida, updated the readme file and checksums file. There were 148 individual samples, so I won’t list them all here.

genefish 23:19 on April 24

Grace’s Notebook: Crab Results progress – GOslim transcriptome plot, macropinocytosis list and notes

Crab transcriptome GOslim

Got the BLAST output from the crab transcriptome to GO slim terms using 042320-bairdi_transcriptome-BLAST-to-GOslim.ipynb

Saved this GOslim text file: GOslim-P-pie.txt

GOslim was made in excel because that was easy for me. Removed “other biological processes” and other non-descriptive GOslim terms.

Macropinocytosis

Macropinocytosis was identified as the only significantly enriched GO term from the 2019 crab dataset comparing the infected and uninfected crabs. (DAVID output: analyses/chart_12EFD.txt)

I made a list of the 12 genes that comprised the macropinocytosis in this Rmd. The list is here.

The genes were all (except one) described in Dictyostelium (slime mold), which is a protist. Hematodinium is also a protist… interesting. I’ve added some notes from papers and other things that I’ve found in the crab paper discussion section.