Last updated: 2024-06-19
Checks: 6 1
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.
Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.
| absolute | relative |
|---|---|
| /Users/alphamanae/Documents/GitHub/locust-comparative-genomics | . |
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 5913857. 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: .RData
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.DS_Store
Ignored: analysis/.Rhistory
Ignored: data/.DS_Store
Ignored: data/.Rhistory
Ignored: data/DEseq2_SAMER_SAMER_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SAMER_SAMER_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SAMER_SGREG_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SAMER_SGREG_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SCANC_SCANC_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SCANC_SCANC_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SCANC_SGREG_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SCANC_SGREG_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SCUBE_SCUBE_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SCUBE_SCUBE_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SCUBE_SGREG_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SCUBE_SGREG_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SGREG_SGREG_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SGREG_SGREG_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SNITE_SGREG_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SNITE_SGREG_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SNITE_SNITE_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SNITE_SNITE_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SPICE_SGREG_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SPICE_SGREG_THORAX_STARnew_features/.DS_Store
Ignored: data/DEseq2_SPICE_SPICE_HEAD_STARnew_features/.DS_Store
Ignored: data/DEseq2_SPICE_SPICE_THORAX_STARnew_features/.DS_Store
Ignored: data/americana/.DS_Store
Ignored: data/americana/deg_counts/.DS_Store
Ignored: data/americana/deg_counts/STAR/.DS_Store
Ignored: data/americana/deg_counts/STAR/americana/.DS_Store
Ignored: data/americana/deg_counts/STAR/americana/featurecounts/.DS_Store
Ignored: data/americana/deg_counts/STAR/gregaria/.DS_Store
Ignored: data/americana/deg_counts/STAR/gregaria/featurecounts/.DS_Store
Ignored: data/americana/deg_counts/STAR_newparams/.DS_Store
Ignored: data/americana/deg_counts/STAR_newparams/americana/.DS_Store
Ignored: data/americana/deg_counts/STAR_newparams/americana/featurecounts/.DS_Store
Ignored: data/americana/deg_counts/STAR_newparams/gregaria/.DS_Store
Ignored: data/americana/deg_counts/STAR_newparams/gregaria/featurecounts/.DS_Store
Ignored: data/cancellata/.DS_Store
Ignored: data/cancellata/deg_counts/.DS_Store
Ignored: data/cancellata/deg_counts/STAR/.DS_Store
Ignored: data/cancellata/deg_counts/STAR_newparams/.DS_Store
Ignored: data/cubense/.DS_Store
Ignored: data/cubense/deg_counts/.DS_Store
Ignored: data/cubense/deg_counts/STAR/.DS_Store
Ignored: data/cubense/deg_counts/STAR/cubense/.DS_Store
Ignored: data/cubense/deg_counts/STAR/cubense/featurecounts/
Ignored: data/cubense/deg_counts/STAR/gregaria/
Ignored: data/cubense/deg_counts/STAR_newparams/.DS_Store
Ignored: data/cubense/deg_counts/STAR_newparams/cubense/.DS_Store
Ignored: data/cubense/deg_counts/STAR_newparams/gregaria/.DS_Store
Ignored: data/gregaria/.DS_Store
Ignored: data/gregaria/deg_counts/.DS_Store
Ignored: data/gregaria/deg_counts/STAR/.DS_Store
Ignored: data/gregaria/deg_counts/STAR_newparams/.DS_Store
Ignored: data/gregaria/deg_counts/STAR_newparams/gregaria/.DS_Store
Ignored: data/gregaria/deg_counts/STAR_newparams/gregaria/featurecounts/.DS_Store
Ignored: data/list/.DS_Store
Ignored: data/metadata/.DS_Store
Ignored: data/nitens/.DS_Store
Ignored: data/nitens/deg_counts/.DS_Store
Ignored: data/nitens/deg_counts/STAR_newparams/.DS_Store
Ignored: data/nitens/deg_counts/STAR_newparams/gregaria/.DS_Store
Ignored: data/nitens/deg_counts/STAR_newparams/gregaria/featurecounts/.DS_Store
Ignored: data/nitens/deg_counts/STAR_newparams/nitens/.DS_Store
Ignored: data/nitens/deg_counts/STAR_newparams/nitens/featurecounts/.DS_Store
Ignored: data/piceifrons/.DS_Store
Ignored: data/piceifrons/deg_counts/.DS_Store
Ignored: data/piceifrons/deg_counts/STAR/.DS_Store
Ignored: data/piceifrons/deg_counts/STAR/gregaria/.DS_Store
Ignored: data/piceifrons/deg_counts/STAR/gregaria/featurecounts/.DS_Store
Ignored: data/piceifrons/deg_counts/STAR/piceifrons/.DS_Store
Ignored: data/piceifrons/deg_counts/STAR/piceifrons/featurecounts/.DS_Store
Ignored: data/piceifrons/deg_counts/STAR_newparams/.DS_Store
Ignored: data/piceifrons/deg_counts/STAR_newparams/gregaria/.DS_Store
Ignored: data/piceifrons/deg_counts/STAR_newparams/piceifrons/.DS_Store
Ignored: figures/
Ignored: tables/
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/1_seq-data.Rmd) and HTML
(docs/1_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 |
|---|---|---|---|---|
| html | 4246194 | Maeva A. TECHER | 2024-01-29 | Build site. |
| html | f701a01 | Maeva A. TECHER | 2024-01-29 | reupdate |
| html | 194fcc7 | Maeva A. TECHER | 2024-01-28 | Build site. |
| Rmd | 811ed09 | Maeva A. TECHER | 2024-01-28 | wflow_publish("analysis/1_seq-data.Rmd") |
| html | 8e9db5e | Maeva A. TECHER | 2024-01-28 | Build site. |
| Rmd | 919028d | Maeva A. TECHER | 2024-01-28 | wflow_publish("analysis/1_seq-data.Rmd") |
| html | 17ad0e9 | Maeva A. TECHER | 2024-01-28 | Build site. |
| Rmd | 895ba26 | Maeva A. TECHER | 2024-01-28 | wflow_publish("analysis/1_seq-data.Rmd") |
| html | 29691fd | Maeva A. TECHER | 2024-01-24 | Build site. |
| html | 1b09cbe | Maeva A. TECHER | 2024-01-24 | remove |
| html | a92c643 | Maeva A. TECHER | 2024-01-24 | Build site. |
| Rmd | 838b0a1 | Maeva A. TECHER | 2024-01-24 | wflow_publish("analysis/1_seq-data.Rmd") |
| Rmd | b195f57 | Maeva A. TECHER | 2024-01-23 | refresh |
| html | b195f57 | Maeva A. TECHER | 2024-01-23 | refresh |
| html | 612efb5 | Maeva A. TECHER | 2023-12-18 | Build site. |
| Rmd | 53877fa | Maeva A. TECHER | 2023-12-18 | add pages |
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("ggplot2")
library("reshape2")
library("plotly")
## 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 Behavioral Plasticity Research Institute (BPRI) successfully generated high-quality, chromosome-length genomes for six species of grasshoppers and locusts within the Schistocerca genus. To achieve this, the genome assembly team utilized a hybrid approach, combining long-read sequencing with HiFi PacBio technology and short-read sequencing with Hi-C. Following the finalization of the genome assemblies, the annotation process was conducted by the Eukaryotic Annotation Pipeline for RefSeq by NCBI.
We will use the RefSeq assemblies that have received accession
numbers starting with “GFC” and the associated annotation .gtf and .gff
files.
Schistocerca
gregaria
Schistocerca
piceifrons
Schistocerca
cancellata
Schistocerca
americana
Schistocerca
serialis cubense
Schistocerca
nitens
Here is an an example of how to download the genome assembly for 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 public sequences with accession on NCBI (DDBJ or other database) or can be generated de novo and yet to be released.
We conducted a search in the NCBI database for RNA SRA associated with the Acrididae, Schistocerca genus and found 199 accessions were available (last update: 25 January 2024). We selected only comparables SRA which were obtained from head and thorax, separately from last instar nymph, females, 3 days after molting.
All RNA-seq data for species including S. americana, S. piceifrons, S. serialis cubense, S. cancellata, S. nitens were originally generated for the BioProject PRJNA633949 at the Song Laboratory of Insect Systematics and Evolution at Texas A&M University and under the same conditions (Foquet et al. 2021).
Using the Run Selector from NCBI, we can easily download a metadata table that allows us to visualize how the accessions are distributed per species. We will use this metadata table for analysis of the RNAseq data from the previously published dataset.
We compile a list of accessions from the “Run” category for
each species and then use SRA-toolkit from NCBI. First, we
create an empty directory named ncbi to download each SRA.
This is where SRA Toolkit will dump the prefetched SRA
files in the .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 sequenced new whole transcriptomes from the desert locust S. gregaria under two different density conditions. While these new transcriptomes originated from specimens reared in similar conditions, we made slight changes in the RNA extraction, library preparation and sequencing processes. Specifically, we used the Illumina Stranded Total RNA with RiboZero depletion kit and sequenced on a NovaSeq SP flow cell at TxGen, aiming for a targeted yield of 40 millions reads per library.
For the sequencing of many transcriptomes at the same time, it is common practice to run pooled libraries across multiples lanes, depending on the required read count. To handle this, I used the following loop to merge several files from 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;We need to create a .csv file containing as much
information as possible for each sample/file name (e.g., Sample_ID,
Species, Sex, RearingCondition). An interactive and searchable table
with this type of information is available below and can even be
downloaded directly.
We can download most of the metadata from NCBI, but as often it is
case, the table is incomplete. I had to look through Biosample,
BioProject description as well as associated articles or grants if
available to complete data especially regarding:
- phase - rearing condition - developmental stage and age - sex -
tissues collected - origin of the RNAseq data from single or pooled
specimens
NB: Throughout our analysis, we will complete this metadata file by adding other stats related to sequencing and mapping.
# Load our SRA metadata table
setwd("/Users/alphamanae/Documents/GitHub/locust-comparative-genomics")
metaseq <- read_table("data/metadata/RNAseq_headthorax_METADATA2024.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"))))We explore quickly how many bulk tissue transcriptomes are available per species, per rearing condition.
# Create a new variable combining Species and RearingCondition
metaseq_modif <- metaseq %>%
mutate(Species_Rearing = paste(Species, RearingCondition, sep = "_"))
# Group and count the data by Species_Rearing and Tissue
plot_data <- metaseq_modif %>%
group_by(Species_Rearing, Tissue, Species) %>%
summarize(Count = n())
# Define colors for each unique Species
species_colors <- c(Schistocerca_gregaria = "orange", Schistocerca_piceifrons = "red",
Schistocerca_cancellata = "purple", Schistocerca_americana = "forestgreen", Schistocerca_serialis_cubense = "yellow",
Schistocerca_nitens = "blue")
# Create the plot with combined variable
gg <- ggplot(plot_data, aes(x = Species_Rearing, y = Tissue, size = Count, fill = Species)) +
geom_point(shape = 21) + labs(x = "Species_Rearing", y = "Tissue", title = "Bubble Chart of Tissue Counts for Species+RearingCondition") +
theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_fill_manual(values = species_colors) + scale_size_continuous(range = c(5,
20))
# Convert the ggplot object to a plotly object for interactivity
interactive_plot <- ggplotly(gg)
# Display the interactive plot
interactive_plot
sessionInfo()FALSE R version 4.3.1 (2023-06-16)
FALSE Platform: x86_64-apple-darwin20 (64-bit)
FALSE Running under: macOS Sonoma 14.4.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_5.1.0 plotly_4.10.4 DT_0.33
FALSE [5] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
FALSE [9] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
FALSE [13] ggplot2_3.5.1 tidyverse_2.0.0 rmdformats_1.0.4 knitr_1.45
FALSE [17] workflowr_1.7.1
FALSE
FALSE loaded via a namespace (and not attached):
FALSE [1] gtable_0.3.5 xfun_0.44 bslib_0.7.0 htmlwidgets_1.6.4
FALSE [5] processx_3.8.4 callr_3.7.6 tzdb_0.4.0 crosstalk_1.2.1
FALSE [9] vctrs_0.6.5 tools_4.3.1 ps_1.7.6 generics_0.1.3
FALSE [13] fansi_1.0.6 highr_0.10 pkgconfig_2.0.3 data.table_1.15.4
FALSE [17] lifecycle_1.0.4 compiler_4.3.1 git2r_0.33.0 munsell_0.5.1
FALSE [21] getPass_0.2-4 httpuv_1.6.15 htmltools_0.5.8.1 sass_0.4.9
FALSE [25] yaml_2.3.8 lazyeval_0.2.2 crayon_1.5.2 later_1.3.2
FALSE [29] pillar_1.9.0 jquerylib_0.1.4 whisker_0.4.1 cachem_1.1.0
FALSE [33] tidyselect_1.2.1 digest_0.6.35 stringi_1.8.4 bookdown_0.39
FALSE [37] labeling_0.4.3 rprojroot_2.0.4 fastmap_1.2.0 grid_4.3.1
FALSE [41] colorspace_2.1-0 cli_3.6.2 magrittr_2.0.3 utf8_1.2.4
FALSE [45] withr_3.0.0 scales_1.3.0 promises_1.3.0 timechange_0.3.0
FALSE [49] rmarkdown_2.27 httr_1.4.7 hms_1.1.3 evaluate_0.23
FALSE [53] viridisLite_0.4.2 rlang_1.1.3 Rcpp_1.0.12 glue_1.7.0
FALSE [57] formatR_1.14 rstudioapi_0.16.0 jsonlite_1.8.8 plyr_1.8.9
FALSE [61] R6_2.5.1 fs_1.6.4