09-Apul-sRNAseq-miRTrace

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

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

---

Inputs:

• Trimmed sRNAseq FastQs generated by 08-Apul-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}/D-Apul/output/09-Apul-sRNAseq-miRTrace'
echo 'export trimmed_reads_dir=${deep_dive_dir}/D-Apul/output/08-Apul-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}/D-Apul/output/09-Apul-sRNAseq-miRTrace
export trimmed_reads_dir=${deep_dive_dir}/D-Apul/output/08-Apul-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/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-140-S1-TP2.flexbar_trim.25bp_1.fastq.gz,sRNA-ACR-140_1
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-140-S1-TP2.flexbar_trim.25bp_2.fastq.gz,sRNA-ACR-140_2
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-145-S1-TP2.flexbar_trim.25bp_1.fastq.gz,sRNA-ACR-145_1
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-145-S1-TP2.flexbar_trim.25bp_2.fastq.gz,sRNA-ACR-145_2
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-150-S1-TP2.flexbar_trim.25bp_1.fastq.gz,sRNA-ACR-150_1
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-150-S1-TP2.flexbar_trim.25bp_2.fastq.gz,sRNA-ACR-150_2
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-173-S1-TP2.flexbar_trim.25bp_1.fastq.gz,sRNA-ACR-173_1
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-173-S1-TP2.flexbar_trim.25bp_2.fastq.gz,sRNA-ACR-173_2
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-178-S1-TP2.flexbar_trim.25bp_1.fastq.gz,sRNA-ACR-178_1
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/trimmed-reads/sRNA-ACR-178-S1-TP2.flexbar_trim.25bp_2.fastq.gz,sRNA-ACR-178_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 10 sample(s).
NOTE: reusing existing output directory, outdated files may be present.

Run complete. Processed 10 sample(s) in 87 s.

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

real    1m28.271s
user    10m24.928s
sys 0m11.388s
/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/09-Apul-sRNAseq-miRTrace
├── [1.6K]  mirtrace.config
├── [302K]  mirtrace-report.html
├── [ 575]  mirtrace-stats-contamination_basic.tsv
├── [ 822]  mirtrace-stats-contamination_detailed.tsv
└── [4.0K]  qc_passed_reads.all.collapsed
    ├── [117M]  sRNA-ACR-140-S1-TP2.flexbar_trim.25bp_1.fasta
    ├── [163M]  sRNA-ACR-140-S1-TP2.flexbar_trim.25bp_2.fasta
    ├── [128M]  sRNA-ACR-145-S1-TP2.flexbar_trim.25bp_1.fasta
    ├── [175M]  sRNA-ACR-145-S1-TP2.flexbar_trim.25bp_2.fasta
    ├── [135M]  sRNA-ACR-150-S1-TP2.flexbar_trim.25bp_1.fasta
    ├── [179M]  sRNA-ACR-150-S1-TP2.flexbar_trim.25bp_2.fasta
    ├── [113M]  sRNA-ACR-173-S1-TP2.flexbar_trim.25bp_1.fasta
    ├── [151M]  sRNA-ACR-173-S1-TP2.flexbar_trim.25bp_2.fasta
    ├── [109M]  sRNA-ACR-178-S1-TP2.flexbar_trim.25bp_1.fasta
    └── [140M]  sRNA-ACR-178-S1-TP2.flexbar_trim.25bp_2.fasta

1 directory, 14 files

4 Results

4.1 Read in table as data frame

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

str(mirtrace.detailed.df)

'data.frame':   7 obs. of  14 variables:
 $ CLADE         : chr  "lophotrochozoa" "lophotrochozoa" "lophotrochozoa" "rodents" ...
 $ FAMILY_ID     : int  1994 1985 1984 351 618 576 576
 $ MIRBASE_IDS   : chr  "cla-miR-1994,cte-miR-1994,lgi-miR-1994a,lgi-miR-1994b" "hru-miR-1985,lgi-miR-1985" "hru-miR-1984,lgi-miR-1984" "mmu-miR-351,rno-miR-351" ...
 $ SEQ           : chr  "TGAGACAGTGTGTCCTCCCT" "TGCCATTTTTATCAGTCACT" "TGCCCTATCCGTCAGGAACT" "TCCCTGAGGAGCCCTTTGAG" ...
 $ sRNA.ACR.140_1: int  0 0 0 0 0 0 0
 $ sRNA.ACR.140_2: int  0 0 0 0 0 0 0
 $ sRNA.ACR.145_1: int  0 0 0 0 0 0 0
 $ sRNA.ACR.145_2: int  0 0 0 0 0 0 0
 $ sRNA.ACR.150_1: int  1 11 10 0 0 1 0
 $ sRNA.ACR.150_2: int  0 0 0 0 0 0 0
 $ sRNA.ACR.173_1: int  3 15 16 0 1 0 1
 $ sRNA.ACR.173_2: int  0 0 0 0 0 0 0
 $ sRNA.ACR.178_1: int  0 0 0 2 0 0 0
 $ sRNA.ACR.178_2: int  0 0 0 0 0 0 0

4.2 Number of samples with matches

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

# 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:  3"

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:  30"

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: 3"

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.ACR.140_1	sRNA.ACR.140_2	sRNA.ACR.145_1	sRNA.ACR.145_2	sRNA.ACR.150_1	sRNA.ACR.150_2	sRNA.ACR.173_1	sRNA.ACR.173_2	sRNA.ACR.178_1	sRNA.ACR.178_2
lophotrochozoa	1994	cla-miR-1994,cte-miR-1994,lgi-miR-1994a,lgi-miR-1994b	TGAGACAGTGTGTCCTCCCT	0	0	0	0	1	0	3	0	0	0
lophotrochozoa	1985	hru-miR-1985,lgi-miR-1985	TGCCATTTTTATCAGTCACT	0	0	0	0	11	0	15	0	0	0
lophotrochozoa	1984	hru-miR-1984,lgi-miR-1984	TGCCCTATCCGTCAGGAACT	0	0	0	0	10	0	16	0	0	0
rodents	351	mmu-miR-351,rno-miR-351	TCCCTGAGGAGCCCTTTGAG	0	0	0	0	0	0	0	0	2	0
primates	618	hsa-miR-618,mml-miR-618,ppy-miR-618,ptr-miR-618	AAACTCTACTTGTCCTTCTG	0	0	0	0	0	0	1	0	0	0
primates	576	hsa-miR-576,mml-miR-576,ppy-miR-576,ptr-miR-576	AAGATGTGGAAAAATTGGAA	0	0	0	0	1	0	0	0	0	0
primates	576	hsa-miR-576	ATTCTAATTTCTCCACGTCT	0	0	0	0	0	0	1	0	0	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.