Last updated: 2023-08-28
Checks: 7 0
Knit directory:
locust-comparative-genomics/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
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.
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.
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.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
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 4044ae4. 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: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.DS_Store
Ignored: analysis/.Rhistory
Ignored: data/.DS_Store
Ignored: data/.Rhistory
Untracked files:
Untracked: data/metadata/
Unstaged changes:
Deleted: data/RNAseq_modified_METADATA2022.txt
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/seq-data.Rmd) and HTML
(docs/seq-data.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 |
|---|---|---|---|---|
| Rmd | 1597fa2 | Maeva A. TECHER | 2023-08-28 | test |
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)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.gzThe 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.
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 --interactiveOnce 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/; doneWe 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.
The details of the sequences are as follow:
Bulk tissue from 2nd generation solitary control
GREG-HATCH-S1-FULL_S20: Solitary hatch-ling 1st instar from Pearl and
Atticus
GREG-S-ICT-9-ALB_S6: Antenna lobes from replicate #9
GREG-S-ICT-9-ANT_S5: Antennae from replicate #9
GREG-S-ICT-9-FAT_S3: Fat body from replicate #9
GREG-S-ICT-9-MHB_S7: Mushroom body from replicate #9
GREG-S-ICT-9-MOP_S4: Maxillary palps from replicate #9
GREG-S-ICT-9-MTG_S2: Metathoracic ganglia from replicate #9
GREG-S-ICT-9-OLB_S8: Optical lobes from replicate #9
GREG-S-ICT-9-WG_S1: Whole gut from replicate #9
Bulk tissue from highly crowded control
GREG-G-CCT-11-ALB_S1: Antenna lobes from replicate #11
GREG-G-CCT-11-ALB-FULL_S17: Antenna lobes from replicate #11
GREG-G-CCT-11-ANT_S13: Antennae from replicate #11
GREG-G-CCT-11-FAT_S11: Fat body from replicate #11
GREG-G-CCT-11-FAT-FULL_S19: Fat body from replicate #11
GREG-G-CCT-11-MHB_S15: Mushroom body from replicate #11
GREG-G-CCT-11-MOP_S12: Maxillary palps from replicate #11
GREG-G-CCT-11-MTG_S10: Metathoracic ganglia from replicate #11
GREG-G-CCT-11-OLB_S16: Optical lobes from replicate #11
GREG-G-CCT-11-OLB-FULL_S18: Optical lobes from replicate #11
GREG-G-CCT-11-WG_S9: Whole gut from replicate #11
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;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.3.1 (2023-06-16)
FALSE Platform: x86_64-apple-darwin20 (64-bit)
FALSE Running under: macOS Ventura 13.5.1
FALSE
FALSE Matrix products: default
FALSE BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
FALSE LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
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 time zone: America/Chicago
FALSE tzcode source: internal
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.2 DT_0.28
FALSE [5] lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.2
FALSE [9] purrr_1.0.2 readr_2.1.4 tidyr_1.3.0 tibble_3.2.1
FALSE [13] ggplot2_3.4.3 tidyverse_2.0.0 rmdformats_1.0.4 knitr_1.43
FALSE [17] workflowr_1.7.1
FALSE
FALSE loaded via a namespace (and not attached):
FALSE [1] gtable_0.3.4 xfun_0.40 bslib_0.5.1 htmlwidgets_1.6.2
FALSE [5] processx_3.8.2 callr_3.7.3 tzdb_0.4.0 crosstalk_1.2.0
FALSE [9] vctrs_0.6.3 tools_4.3.1 ps_1.7.5 generics_0.1.3
FALSE [13] fansi_1.0.4 pkgconfig_2.0.3 data.table_1.14.8 lifecycle_1.0.3
FALSE [17] compiler_4.3.1 git2r_0.32.0 munsell_0.5.0 getPass_0.2-2
FALSE [21] httpuv_1.6.11 htmltools_0.5.6 sass_0.4.7 yaml_2.3.7
FALSE [25] lazyeval_0.2.2 crayon_1.5.2 later_1.3.1 pillar_1.9.0
FALSE [29] jquerylib_0.1.4 whisker_0.4.1 ellipsis_0.3.2 cachem_1.0.8
FALSE [33] tidyselect_1.2.0 digest_0.6.33 stringi_1.7.12 bookdown_0.35
FALSE [37] rprojroot_2.0.3 fastmap_1.1.1 grid_4.3.1 colorspace_2.1-0
FALSE [41] cli_3.6.1 magrittr_2.0.3 utf8_1.2.3 withr_2.5.0
FALSE [45] scales_1.2.1 promises_1.2.1 timechange_0.2.0 rmarkdown_2.24
FALSE [49] httr_1.4.7 hms_1.1.3 evaluate_0.21 viridisLite_0.4.2
FALSE [53] rlang_1.1.1 Rcpp_1.0.11 glue_1.6.2 formatR_1.14
FALSE [57] rstudioapi_0.15.0 jsonlite_1.8.7 R6_2.5.1 plyr_1.8.8
FALSE [61] fs_1.6.3