09-Peve-sRNAseq-miRTrace

Use miRTrace (Kang et al. 2018) to identify taxonomic origins of miRNA sequencing data.

NOTE: This requires you to have previously run 06-Peve-sRNAseq-trimming.Rmd, as the code relies on the trimmed reads output from that code.

---

Inputs:

• Trimmed sRNAseq FastQs generated by 06-Peve-sRNAseq-trimming.Rmd

  • Filenames formatted: *flexbar_trim.25bp*.gz

Outputs:

• mirtrace.config: A miRTrace config file. A comma-separated file with this layout (one FastQ per line): /path/to/fastq,custom_sample_name

• “Collapsed” (i.e. unique sequences only) FastA for each corresponding input FastQ.

• mirtrace-report.html: HTML-formatted report generated by miRTrace.

• mirtrace-stats-contamination_basic.tsv: Tab-delimited report with counts of sequences from each collapsed FastAs having matches to known miRNAs within each of the miRTrace Clades.

• mirtrace-stats-contamination_detailed.csv: Tab-delimited report of only Clades with which sequences were matched, along with the corresponding miRNA families in each clade, and the sequence counts.

1 Create a Bash variables file

This allows usage of Bash variables across R Markdown chunks.

{
echo "#### Assign Variables ####"
echo ""

echo "# Data directories"
echo 'export deep_dive_dir=/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive'
echo 'export output_dir_top=${deep_dive_dir}/E-Peve/output/09-Peve-sRNAseq-miRTrace'
echo 'export trimmed_reads_dir=${deep_dive_dir}/E-Peve/output/06-Peve-sRNAseq-trimming/trimmed-reads'
echo ""

echo "# Paths to programs"
echo 'export mirtrace=/home/sam/programs/mambaforge/envs/miRTrace_env/bin/mirtrace'
echo ""

echo "# Set number of CPUs to use"
echo 'export threads=40'
echo ""

echo "export fastq_pattern='*flexbar_trim.25bp*.gz'"

echo "# Programs associative array"
echo "declare -A programs_array"
echo "programs_array=("
echo '[mirtrace]="${mirtrace}"'
echo ")"
} > .bashvars

cat .bashvars

#### Assign Variables ####

# Data directories
export deep_dive_dir=/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive
export output_dir_top=${deep_dive_dir}/E-Peve/output/09-Peve-sRNAseq-miRTrace
export trimmed_reads_dir=${deep_dive_dir}/E-Peve/output/06-Peve-sRNAseq-trimming/trimmed-reads

# Paths to programs
export mirtrace=/home/sam/programs/mambaforge/envs/miRTrace_env/bin/mirtrace

# Set number of CPUs to use
export threads=40

export fastq_pattern='*flexbar_trim.25bp*.gz'
# Programs associative array
declare -A programs_array
programs_array=(
[mirtrace]="${mirtrace}"
)

2 Create miRTrace config file

# Load bash variables into memory
source .bashvars

# Declare array
fastq_array=()

# Populate array
fastq_array=(${trimmed_reads_dir}/${fastq_pattern})

# Loop through read pairs
# Increment by 2 to process next pair of FastQ files
if [ -f "${output_dir_top}/mirtrace.config" ]; then
  echo "mirtrace.config already exists. Nothing to do."
  
else

  for (( i=0; i<${#fastq_array[@]} ; i+=2 ))
  do
    # Use first three parts of filename to create short sample name
    R1_name=$(echo "${fastq_array[i]##*/}" | awk -F "-" '{print $1"-"$2"-"$3}')
    R2_name=$(echo "${fastq_array[i+1]##*/}" | awk -F "-" '{print $1"-"$2"-"$3}')
    echo "${fastq_array[i]},${R1_name}_1"
    echo "${fastq_array[i+1]},${R2_name}_2"
  done >> "${output_dir_top}/mirtrace.config"

fi

cat "${output_dir_top}/mirtrace.config"

/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/06-Peve-sRNAseq-trimming/trimmed-reads/sRNA-POR-73-S1-TP2.flexbar_trim.25bp_1.fastq.gz,sRNA-POR-73_1
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/06-Peve-sRNAseq-trimming/trimmed-reads/sRNA-POR-73-S1-TP2.flexbar_trim.25bp_2.fastq.gz,sRNA-POR-73_2
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/06-Peve-sRNAseq-trimming/trimmed-reads/sRNA-POR-79-S1-TP2.flexbar_trim.25bp_1.fastq.gz,sRNA-POR-79_1
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/06-Peve-sRNAseq-trimming/trimmed-reads/sRNA-POR-79-S1-TP2.flexbar_trim.25bp_2.fastq.gz,sRNA-POR-79_2
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/06-Peve-sRNAseq-trimming/trimmed-reads/sRNA-POR-82-S1-TP2.flexbar_trim.25bp_1.fastq.gz,sRNA-POR-82_1
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/06-Peve-sRNAseq-trimming/trimmed-reads/sRNA-POR-82-S1-TP2.flexbar_trim.25bp_2.fastq.gz,sRNA-POR-82_2

3 Run miRTrace

# Load bash variables into memory
source .bashvars

time \
${programs_array[mirtrace]} trace \
--config ${output_dir_top}/mirtrace.config \
--write-fasta \
--num-threads ${threads} \
--output-dir ${output_dir_top} \
--force

tree -h ${output_dir_top}

miRTrace version 1.0.1 starting. Processing 6 sample(s).
NOTE: reusing existing output directory, outdated files may be present.

Run complete. Processed 6 sample(s) in 38 s.

Reports written to /home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/09-Peve-sRNAseq-miRTrace/
For information about citing our paper, run miRTrace in mode "cite".

real    0m39.313s
user    2m46.272s
sys 0m4.148s
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/E-Peve/output/09-Peve-sRNAseq-miRTrace
├── [ 972]  mirtrace.config
├── [299K]  mirtrace-report.html
├── [ 395]  mirtrace-stats-contamination_basic.tsv
├── [ 413]  mirtrace-stats-contamination_detailed.tsv
└── [4.0K]  qc_passed_reads.all.collapsed
    ├── [ 63M]  sRNA-POR-73-S1-TP2.flexbar_trim.25bp_1.fasta
    ├── [ 74M]  sRNA-POR-73-S1-TP2.flexbar_trim.25bp_2.fasta
    ├── [ 94M]  sRNA-POR-79-S1-TP2.flexbar_trim.25bp_1.fasta
    ├── [109M]  sRNA-POR-79-S1-TP2.flexbar_trim.25bp_2.fasta
    ├── [112M]  sRNA-POR-82-S1-TP2.flexbar_trim.25bp_1.fasta
    └── [128M]  sRNA-POR-82-S1-TP2.flexbar_trim.25bp_2.fasta

1 directory, 10 files

4 Results

4.1 Read in table as data frame

mirtrace.detailed.df <- read.csv("../output/09-Peve-sRNAseq-miRTrace/mirtrace-stats-contamination_detailed.tsv", sep = "\t", header = TRUE)

str(mirtrace.detailed.df)

'data.frame':   1 obs. of  10 variables:
 $ CLADE        : chr "insects"
 $ FAMILY_ID    : int 14
 $ MIRBASE_IDS  : chr "aae-miR-14,aga-miR-14,ame-miR-14,api-miR-14,bmo-miR-14,cqu-miR-14,dan-miR-14,der-miR-14,dgr-miR-14,dme-miR-14,d"| __truncated__
 $ SEQ          : chr "TCAGTCTTTTTCTCTCTCCT"
 $ sRNA.POR.73_1: int 0
 $ sRNA.POR.73_2: int 0
 $ sRNA.POR.79_1: int 0
 $ sRNA.POR.79_2: int 0
 $ sRNA.POR.82_1: int 1
 $ sRNA.POR.82_2: int 0

4.2 Number of samples with matches

IMPORTANT: Change starts_with() to match prefix of input sRNAseq reads!

# Select columns corresponding to sample names
sample_columns <- mirtrace.detailed.df %>%
  select(starts_with("sRNA.POR."))

# Calculate the sum for each column
sample_sums <- colSums(sample_columns)

# Count the number of columns with a sum greater than 0
samples_with_sum_gt_0 <- sum(sample_sums > 0)

paste("Number of samples with matches: ", samples_with_sum_gt_0)

[1] "Number of samples with matches:  1"

4.3 Percentage of samples with matches

# Total number of samples (columns)
total_samples <- ncol(sample_columns)

# Percentage of samples with sums greater than 0
percentage_samples_gt_0 <- (samples_with_sum_gt_0 / total_samples) * 100

paste("Percentage of samples with matches: ", percentage_samples_gt_0)

[1] "Percentage of samples with matches:  16.6666666666667"

4.4 Number of clades with matches

unique_clade_count <- mirtrace.detailed.df %>%
  distinct(CLADE) %>%    # Get unique entries in CLADE column
  count()               # Count the number of unique entries



paste("Number of clades with matches:", unique_clade_count)

[1] "Number of clades with matches: 1"

4.5 miRTrace table

To make them easier to see, counts > 0 are highlighted in green.

mirtrace.detailed.df %>%
  mutate(
    across(
      starts_with("sRNA"),
      ~cell_spec(
        .,
        background = ifelse(
          . > 0,
          "lightgreen",
          "white"
          )
        )
      )
    ) %>%
  kable(escape = F, caption = "Clades identified as having sRNAseq matches.") %>%
  kable_styling("striped") %>% 
  scroll_box(width = "100%", height = "500px")

Table 4.1: Clades identified as having sRNAseq matches.

CLADE	FAMILY_ID	MIRBASE_IDS	SEQ	sRNA.POR.73_1	sRNA.POR.73_2	sRNA.POR.79_1	sRNA.POR.79_2	sRNA.POR.82_1	sRNA.POR.82_2
insects	14	aae-miR-14,aga-miR-14,ame-miR-14,api-miR-14,bmo-miR-14,cqu-miR-14,dan-miR-14,der-miR-14,dgr-miR-14,dme-miR-14,dmo-miR-14,dpe-miR-14,dps-miR-14,dse-miR-14,dsi-miR-14,dvi-miR-14,dwi-miR-14,dya-miR-14,hme-miR-14,mse-miR-14,ngi-miR-14,nvi-miR-14,tca-miR-14	TCAGTCTTTTTCTCTCTCCT	0	0	0	0	1	0

---

Citations

Kang, Wenjing, Yrin Eldfjell, Bastian Fromm, Xavier Estivill, Inna Biryukova, and Marc R. Friedländer. 2018. “miRTrace Reveals the Organismal Origins of microRNA Sequencing Data.” Genome Biology 19 (1). https://doi.org/10.1186/s13059-018-1588-9.