Last updated: 2022-10-31

Checks: 7 0

Knit directory: locust-phase-transition-RNAseq/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20221025)

The command set.seed(20221025) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 2ccb4bc

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 2ccb4bc. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    analysis/.DS_Store
    Ignored:    data/.DS_Store
    Ignored:    data/americana/.DS_Store
    Ignored:    data/americana/STAR_counts_4thcol/.DS_Store
    Ignored:    data/cancellata/.DS_Store
    Ignored:    data/cancellata/STAR_counts_4thcol/.DS_Store
    Ignored:    data/cubense/.DS_Store
    Ignored:    data/cubense/STAR_counts_4thcol/.DS_Store
    Ignored:    data/gregaria/.DS_Store
    Ignored:    data/gregaria/STAR_counts_4thcol/.DS_Store
    Ignored:    data/gregaria/list/.DS_Store
    Ignored:    data/metadata/.DS_Store
    Ignored:    data/nitens/.DS_Store
    Ignored:    data/nitens/STAR_counts_4thcol/.DS_Store
    Ignored:    data/piceifrons/.DS_Store
    Ignored:    data/piceifrons/STAR_counts_4thcol/.DS_Store
    Ignored:    data/piceifrons/list/.DS_Store

Unstaged changes:
    Deleted:    analysis/gene-quant.Rmd
    Deleted:    analysis/own-refseq.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/seqdata-prep.Rmd) and HTML (docs/seqdata-prep.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
html	8481823	MaevaTecher	2022-10-31	Build site.
html	ea70380	MaevaTecher	2022-10-31	Build site.
Rmd	de64b6b	MaevaTecher	2022-10-31	wflow_publish(c("analysis/_site.yml", "analysis/map-refseq.Rmd",
html	8b34f19	MaevaTecher	2022-10-30	Build site.
html	bd04bb5	MaevaTecher	2022-10-30	Build site.
html	61f4b75	MaevaTecher	2022-10-30	Build site.
Rmd	f65c630	MaevaTecher	2022-10-30	wflow_publish(c("analysis/_site.yml", "analysis/gene-quant.Rmd",

Load R libraries (install first from CRAN or Bioconductor)

library("knitr")
library("rmdformats")
library("tidyverse")
library("DT")  # for making interactive search table
library("plotly") # for interactive plots
library("ggthemes") # for theme_calc
library("reshape2")

## Global options
options(max.print="10000")
knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE,
    cache = FALSE,
    comment = FALSE,
    prompt = FALSE,
    tidy = TRUE
)
opts_knit$set(width=75)

Schistocerca reference genome

The BPRI generated high-quality and chromosome length genome for six species of grasshoppers and locusts Schistocerca. For this, the genome assembly team used a combination of long-reads sequencing with HiFi PacBio and short-reads with Hi-C. Following the finalization of the genomes assemblies, the annotation was conducted by the Eukaryotic Annotation Pipeline for RefSeq by NCBI.

We will use the assemblies that received accession numbers for RefSeq and associated with .gtf and .gff files.

Below is an example for downloading Schistocerca piceifrons from NCBI RefSeq:

## Downloading the genome sequence, primary assembly fasta file with RefSeq contig accessions
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/021/461/385/GCF_021461385.2_iqSchPice1.1/GCF_021461385.2_iqSchPice1.1_genomic.fna.gz

## Downloading the comprehensive annotation file in .gtf format
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/021/461/385/GCF_021461385.2_iqSchPice1.1/GCF_021461385.2_iqSchPice1.1_genomic.gtf.gz

Obtaining RNA sequence reads

The reads that can be mapped on these sequences can be either from sequences with accession on NCBI, DDBJ or yet to be release after being generated de novo.

Downloading NCBI Transcriptome SRAs

We searched through the NCBI database for RNA SRA associated with Schistocerca in general and found 160 accessions were available. All RNA-seq data for S. americana, S. piceifrons, S. cubense, S. cancellata, S. nitens were generated originally in the Song Lab at Texas A&M University and were produced following the same standard protocol.

Using the Run Selector from NCBI, we can easily download a metadata table which we can use to visualize how the accessions are distributed per species. We will use this metadata table for analysis in which we include pre-generated SRA.

We collect a list of accessions from Run for each species and then use SRA-toolkit from NCBI. First, we make an empty directory named ncbi to download each SRA. This is where SRA Toolkit will dump the prefetched SRA files in .sra format.

ml purge
ml GCC/10.2.0 OpenMPI/4.0.5 SRA-Toolkit/2.10.9
vdb-config --interactive

Once in the vdb-config interactive mode, select cache, choose, then use [ .. ], to enter /home/USERNAME/PATH/ncbi one directory at a time

prefetch --option-file SraAccList.txt
cat SraAccList.txt | xargs fasterq-dump --split-3 --outdir "/your-directory/for-fastq"

Clean-up the ncbi directory and move the fastq.gz file (rename if wanted).

another option

for x in *.sra ; do fasterq-dump --split-files $x ; mv *.fastq ../../paired_end_piceifrons/; done

De novo sequencing

We generated new whole transcriptomes from two density conditions for the desert locust (Acrididae: Schistocerca gregaria). Here, we analyze 20 transcriptomes for a pilot project using Illumina Stranded Total RNA with RiboZero depletion and sequenced on a NovaSeq SP flow cell at TxGen for a targeted amount of 40 millions reads per library.

Merging several runs into a single fastq.gz file

For transcriptome that have yet to be released, it is common for libraries to be run on multiples lanes depending on the amount of reads needed. I used the following loop to merge several files with different Lanes:

#!/bin/bash
# add the print $3 for TxGen reads as it is for samples S
for i in $(ls -1 *R1*.gz | awk -F '_' '{print $1"_"$2"_"$3}' | sort | uniq)
do echo $i
echo "Merging R1 ${i}"
cat "$i"_L00*_R1_001.fastq.gz > "$i"_MERGE_R1_001.fastq.gz
echo "Merging R2 ${i}"
cat "$i"_L00*_R2_001.fastq.gz > "$i"_MERGE_R2_001.fastq.gz
done;

Preparing the metadata file

Make a .csv file with as much information as possible per sample/file name (e.g., Sample_ID, Species, Sex, RearingCondition). An interactive and searchable table is found below and even be downloaded directly.

NB: Throughout our analysis, we will complete this metadata file by adding other stats related to sequencing and mapping.

# Load our SRA metadata table
metaseq <- read_table("data/metadata/RNAseq_modified_METADATA2022.txt", col_names = TRUE,
    guess_max = 5000)

## Create an interactive search table
metaseq %>%
    datatable(extensions = "Buttons", options = list(dom = "Blfrtip", buttons = c("copy",
        "csv", "excel"), lengthMenu = list(c(10, 20, 50, 100, 200, -1), c(10, 20,
        50, 100, 200, "All"))))

sessionInfo()

FALSE R version 4.2.1 (2022-06-23)
FALSE Platform: x86_64-apple-darwin17.0 (64-bit)
FALSE Running under: macOS Big Sur ... 10.16
FALSE 
FALSE Matrix products: default
FALSE BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
FALSE LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
FALSE 
FALSE locale:
FALSE [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
FALSE 
FALSE attached base packages:
FALSE [1] stats     graphics  grDevices utils     datasets  methods   base     
FALSE 
FALSE other attached packages:
FALSE  [1] reshape2_1.4.4   ggthemes_4.2.4   plotly_4.10.0    DT_0.26         
FALSE  [5] forcats_0.5.2    stringr_1.4.1    dplyr_1.0.10     purrr_0.3.5     
FALSE  [9] readr_2.1.3      tidyr_1.2.1      tibble_3.1.8     ggplot2_3.3.6   
FALSE [13] tidyverse_1.3.2  rmdformats_1.0.4 knitr_1.40       workflowr_1.7.0 
FALSE 
FALSE loaded via a namespace (and not attached):
FALSE  [1] fs_1.5.2            lubridate_1.8.0     httr_1.4.4         
FALSE  [4] rprojroot_2.0.3     tools_4.2.1         backports_1.4.1    
FALSE  [7] bslib_0.4.0         utf8_1.2.2          R6_2.5.1           
FALSE [10] DBI_1.1.3           lazyeval_0.2.2      colorspace_2.0-3   
FALSE [13] withr_2.5.0         tidyselect_1.2.0    processx_3.7.0     
FALSE [16] compiler_4.2.1      git2r_0.30.1        cli_3.4.1          
FALSE [19] rvest_1.0.3         formatR_1.12        xml2_1.3.3         
FALSE [22] bookdown_0.29       sass_0.4.2          scales_1.2.1       
FALSE [25] callr_3.7.2         digest_0.6.30       rmarkdown_2.17     
FALSE [28] pkgconfig_2.0.3     htmltools_0.5.3     dbplyr_2.2.1       
FALSE [31] fastmap_1.1.0       htmlwidgets_1.5.4   rlang_1.0.6        
FALSE [34] readxl_1.4.1        rstudioapi_0.14     jquerylib_0.1.4    
FALSE [37] generics_0.1.3      jsonlite_1.8.3      crosstalk_1.2.0    
FALSE [40] googlesheets4_1.0.1 magrittr_2.0.3      Rcpp_1.0.9         
FALSE [43] munsell_0.5.0       fansi_1.0.3         lifecycle_1.0.3    
FALSE [46] stringi_1.7.8       whisker_0.4         yaml_2.3.6         
FALSE [49] plyr_1.8.7          grid_4.2.1          promises_1.2.0.1   
FALSE [52] crayon_1.5.2        haven_2.5.1         hms_1.1.2          
FALSE [55] ps_1.7.1            pillar_1.8.1        reprex_2.0.2       
FALSE [58] glue_1.6.2          evaluate_0.17       getPass_0.2-2      
FALSE [61] data.table_1.14.4   modelr_0.1.9        vctrs_0.5.0        
FALSE [64] tzdb_0.3.0          httpuv_1.6.6        cellranger_1.1.0   
FALSE [67] gtable_0.3.1        assertthat_0.2.1    cachem_1.0.6       
FALSE [70] xfun_0.34           broom_1.0.1         later_1.3.0        
FALSE [73] googledrive_2.0.0   viridisLite_0.4.1   gargle_1.2.1       
FALSE [76] ellipsis_0.3.2

Downloading or obtaining sequencing data

Maeva Techer

2022-10-31