13.1-Apul-sRNAseq-ShortStack-R1-reads

Use ShortStack (Axtell 2013; Shahid and Axtell 2014; Johnson et al. 2016)to perform alignment of sRNAseq data and annotation of sRNA-producing genes.

The A.millepora genome will be used as the reference genome for A.pulchra, as A.pulchra does not currently have a sequenced genome and A.millepora had highest alignment rates for standard RNAseq data compared to other published genomes tested.

---

Inputs:

• Requires trimmed sRNAseq files generated by 08.1-Dapul-sRNAseq-trimming-R1-only.Rmd

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

• A.millepora genome FastA. See 12-Apul-sRNAseq-MirMachine.Rmd for download info if needed.

Outputs:

• See ShortStack outputs documentation for full list and detailed descriptions.

Software requirements:

• Utilizes a ShortStack Conda/Mamba environment, per the installation instructions.

Replace with name of your ShortStack environment and the path to the corresponding conda installation (find this after you’ve activated the environment).

E.g.

# Activate environment
conda activate ShortStack4_env

# Find conda path
which conda

---

1 Set R variables

shortstack_conda_env_name <- c("ShortStack-4.0.3_env")
shortstack_cond_path <- c("/home/sam/programs/mambaforge/condabin/conda")

2 Create a Bash variables file

This allows usage of Bash variables across R Markdown chunks.

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

echo "# Trimmed FastQ naming pattern"
echo "export trimmed_fastqs_pattern='*fastp-R1-31bp-auto_adapters-polyG.fq.gz'"

echo "# Data directories"
echo 'export deep_dive_dir=/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive'
echo 'export deep_dive_data_dir="${deep_dive_dir}/data"'
echo 'export output_dir_top=${deep_dive_dir}/D-Apul/output/13.1-Apul-sRNAseq-ShortStack-R1-reads'
echo 'export trimmed_fastqs_dir="${deep_dive_dir}/D-Apul/output/08.1-Dapul-sRNAseq-trimming-R1-only/trimmed-reads"'
echo ""

echo "# Input/Output files"
echo 'export genome_fasta_dir=${deep_dive_dir}/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1'
echo 'export genome_fasta_name="GCF_013753865.1_Amil_v2.1_genomic.fna"'
echo 'export shortstack_genome_fasta_name="GCF_013753865.1_Amil_v2.1_genomic.fa"'
echo 'export mirbase_mature_fasta=mature.fa'
echo 'export mirbase_mature_fasta_version=mirbase-mature-v22.1.fa'
echo 'export genome_fasta="${genome_fasta_dir}/${shortstack_genome_fasta_name}"'
echo ""

echo "# External data URLs"
echo 'export mirbase_fasta_url="https://mirbase.org/download_version_files/22.1/"'
echo ""

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

echo "# Initialize arrays"
echo 'export trimmed_fastqs_array=()'


} > .bashvars

cat .bashvars

#### Assign Variables ####

# Trimmed FastQ naming pattern
export trimmed_fastqs_pattern='*fastp-R1-31bp-auto_adapters-polyG.fq.gz'
# Data directories
export deep_dive_dir=/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive
export deep_dive_data_dir="${deep_dive_dir}/data"
export output_dir_top=${deep_dive_dir}/D-Apul/output/13.1-Apul-sRNAseq-ShortStack-R1-reads
export trimmed_fastqs_dir="${deep_dive_dir}/D-Apul/output/08.1-Dapul-sRNAseq-trimming-R1-only/trimmed-reads"

# Input/Output files
export genome_fasta_dir=${deep_dive_dir}/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1
export genome_fasta_name="GCF_013753865.1_Amil_v2.1_genomic.fna"
export shortstack_genome_fasta_name="GCF_013753865.1_Amil_v2.1_genomic.fa"
export mirbase_mature_fasta=mature.fa
export mirbase_mature_fasta_version=mirbase-mature-v22.1.fa
export genome_fasta="${genome_fasta_dir}/${shortstack_genome_fasta_name}"

# External data URLs
export mirbase_fasta_url="https://mirbase.org/download_version_files/22.1/"

# Set number of CPUs to use
export threads=40

# Initialize arrays
export trimmed_fastqs_array=()

3 Load ShortStack conda environment

If this is successful, the first line of output should show that the Python being used is the one in your [ShortStack](https://github.com/MikeAxtell/ShortStack conda environment path.

E.g.

python: /home/sam/programs/mambaforge/envs/mirmachine_env/bin/python

use_condaenv(condaenv = shortstack_conda_env_name, conda = shortstack_cond_path)

# Check successful env loading
py_config()

python:         /home/sam/programs/mambaforge/envs/ShortStack-4.0.3_env/bin/python
libpython:      /home/sam/programs/mambaforge/envs/ShortStack-4.0.3_env/lib/libpython3.10.so
pythonhome:     /home/sam/programs/mambaforge/envs/ShortStack-4.0.3_env:/home/sam/programs/mambaforge/envs/ShortStack-4.0.3_env
version:        3.10.13 | packaged by conda-forge | (main, Dec 23 2023, 15:36:39) [GCC 12.3.0]
numpy:          /home/sam/programs/mambaforge/envs/ShortStack-4.0.3_env/lib/python3.10/site-packages/numpy
numpy_version:  1.26.4

NOTE: Python version was forced by use_python() function

4 Download miRBase mature miRNA FastA

# Load bash variables into memory
source .bashvars

wget \
--directory-prefix ${deep_dive_data_dir} \
--recursive \
--no-check-certificate \
--continue \
--no-host-directories \
--no-directories \
--no-parent \
--quiet \
--execute robots=off \
 ${mirbase_fasta_url}/${mirbase_mature_fasta}

# Rename to indicate miRBase FastA version
mv ${deep_dive_data_dir}/${mirbase_mature_fasta} ${deep_dive_data_dir}/${mirbase_mature_fasta_version}

ls -lh "${deep_dive_data_dir}"

total 13M
drwxr-xr-x 2 sam sam 4.0K Nov  7 14:36 blast_dbs
-rw-r--r-- 1 sam sam 3.7M Feb 14 16:46 mirbase-mature-v22.1.fa
-rw-r--r-- 1 sam sam 3.7M Dec  4 11:05 mirbase-mature-v22.1-no_spaces.fa
-rw-r--r-- 1 sam sam 3.7M Nov 17 07:59 mirbase-mature-v22.1-no_U.fa
-rw-r--r-- 1 sam sam 726K Nov  7 14:36 mirgene-mature-all-v2.1.fa
-rw-r--r-- 1 sam sam 726K Nov 17 07:59 mirgene-mature-all-v2.1-no_U.fa

5 Run ShortStack

5.1 Modify genome filename for ShortStack compatability

# Load bash variables into memory
source .bashvars

# Check for FastA file first
# Then create rename file if doesn't exist
if [ -f "${genome_fasta_dir}/${shortstack_genome_fasta_name}" ]; then
  echo "${genome_fasta_dir}/${shortstack_genome_fasta_name} already exists. Nothing to do."
  echo ""
else

  # Copy genome FastA to ShortStack-compatible filename (ending with .fa)
  cp ${genome_fasta_dir}/${genome_fasta_name} ${genome_fasta_dir}/${shortstack_genome_fasta_name}
fi

# Confirm
ls -lh ${genome_fasta_dir}/${shortstack_genome_fasta_name}

/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/GCF_013753865.1_Amil_v2.1_genomic.fa already exists. Nothing to do.

-rw-r--r-- 1 sam sam 460M Nov  6 12:40 /home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/data/Amil/ncbi_dataset/data/GCF_013753865.1/GCF_013753865.1_Amil_v2.1_genomic.fa

5.2 Excecute ShortStack command

Uses the --dn_mirna option to identify miRNAs in the genome, without relying on the --known_miRNAs.

This part of the code redirects the output of time to the end of shortstack.log file.

• ; } \ 2>> ${output_dir_top}/shortstack.log

# Load bash variables into memory
source .bashvars

# Make output directory, if it doesn't exist
mkdir --parents "${output_dir_top}"

# Create array of trimmed FastQs
trimmed_fastqs_array=(${trimmed_fastqs_dir}/${trimmed_fastqs_pattern})


# Pass array contents to new variable as space-delimited list
trimmed_fastqs_list=$(echo "${trimmed_fastqs_array[*]}")


###### Run ShortStack ######
{ time \
ShortStack \
--genomefile "${genome_fasta}" \
--readfile ${trimmed_fastqs_list} \
--known_miRNAs ${deep_dive_data_dir}/${mirbase_mature_fasta_version} \
--dn_mirna \
--threads ${threads} \
--outdir ${output_dir_top}/ShortStack_out \
&> ${output_dir_top}/shortstack.log ; } \
2>> ${output_dir_top}/shortstack.log

5.3 Check runtime

# Load bash variables into memory
source .bashvars

tail -n 3 ${output_dir_top}/shortstack.log \
| grep "real" \
| awk '{print "ShortStack runtime:" "\t" $2}'

ShortStack runtime: 52m36.154s

6 Results

6.1 ShortStack synopsis

# Load bash variables into memory
source .bashvars

tail -n 25 ${output_dir_top}/shortstack.log

Writing final files

Found a total of 27 MIRNA loci


Non-MIRNA loci by DicerCall:
N 19203
22 37
23 33
21 9
24 6

Creating visualizations of microRNA loci with strucVis
<<< WARNING >>>
Do not rely on these results alone to annotate new MIRNA loci!
The false positive rate for de novo MIRNA identification is low, but NOT ZERO
Insepct each mirna locus, especially the strucVis output, and see
https://doi.org/10.1105/tpc.17.00851 , https://doi.org/10.1093/nar/gky1141

Wed 14 Feb 2024 16:18:45 -0800 PST
Run Completed!

real    52m36.154s
user    934m28.985s
sys 293m32.838s

ShortStack identified 27 miRNAs.

6.2 Inspect Results.txt

# Load bash variables into memory
source .bashvars

head ${output_dir_top}/ShortStack_out/Results.txt

echo ""
echo "----------------------------------------------------------"
echo ""

echo "Nummber of potential loci:"
awk '(NR>1)' ${output_dir_top}/ShortStack_out/Results.txt | wc -l

Locus   Name    Chrom   Start   End Length  Reads   DistinctSequences   FracTop Strand  MajorRNA    MajorRNAReads   Short   Long    21  22  23  24  DicerCall   MIRNA   known_miRNAs
NC_058066.1:152483-152910   Cluster_1   NC_058066.1 152483  152910  428 142 33  0.04929577464788732 -   UAAGUACUUUAUCAACUAACUCUAGGCA    77  1   130 0   2   0   9   N   N   NA
NC_058066.1:161064-161674   Cluster_2   NC_058066.1 161064  161674  611 517 229 0.29593810444874274 .   UUUUAGCCUAGUGCGGGUUUCCAGACGU    42  32  446 13  5   2   19  N   N   NA
NC_058066.1:172073-172496   Cluster_3   NC_058066.1 172073  172496  424 103 37  0.0970873786407767  -   GCGAUUAUUAACGGCUGGAACGACAGGCG   14  3   82  0   1   0   17  N   N   NA
NC_058066.1:203241-203651   Cluster_4   NC_058066.1 203241  203651  411 107 41  0.6448598130841121  .   UUCUGACUCUAUUAGCAACGAAGACUUU    34  0   106 0   1   0   0   N   N   NA
NC_058066.1:204533-205150   Cluster_5   NC_058066.1 204533  205150  618 326 173 0.7638036809815951  .   UCCCAACACGUCUAGACUGUACAAUUUCU   30  2   312 2   2   5   3   N   N   NA
NC_058066.1:205742-206966   Cluster_6   NC_058066.1 205742  206966  1225    1941    432 0.3286965481710459  .   CAAAAGAGCGGACAAAAUAGUCGACAGAUU  800 9   1886    11  6   9   20  N   N   NA
NC_058066.1:210841-211344   Cluster_7   NC_058066.1 210841  211344  504 1236    347 0.7354368932038835  .   UAAUACUUGUAGUGAAGGUUCAAUCUCGA   90  13  1113    6   11  21  72  N   N   NA
NC_058066.1:349655-351297   Cluster_8   NC_058066.1 349655  351297  1643    3402    1236    0.8080540858318636  +   UCAGCUUGGAAAUGACAGCUUUUGACGU    286 49  3234    10  24  24  61  N   N   NA
NC_058066.1:351491-353439   Cluster_9   NC_058066.1 351491  353439  1949    8987    1695    0.4189384666740848  .   UUUCAAAUCAAAGAUCUUCGCAACGAUGA   766 112 8519    31  47  130 148 N   N   NA

----------------------------------------------------------

Nummber of potential loci:
19315

Column 20 of the Results.txt file identifies if a cluster is a miRNA or not (Y or N).

# Load bash variables into memory
source .bashvars

echo "Number of loci characterized as miRNA:"
awk '$20=="Y" {print $0}' ${output_dir_top}/ShortStack_out/Results.txt \
| wc -l
echo ""

echo "----------------------------------------------------------"

echo ""
echo "Number of loci _not_ characterized as miRNA:"
awk '$20=="N" {print $0}' ${output_dir_top}/ShortStack_out/Results.txt \
| wc -l

Number of loci characterized as miRNA:
27

----------------------------------------------------------

Number of loci _not_ characterized as miRNA:
19288

Column 21 of the Results.txt file identifies if a cluster aligned to a known miRNA (miRBase) or not (Y or NA).

Since there are no miRNAs, the following code will not print any output.

The echo command after the awk command is simply there to prove that the chunk executed.

# Load bash variables into memory
source .bashvars

echo "Number of loci matching miRBase miRNAs:"
awk '$21!="NA" {print $0}' ${output_dir_top}/ShortStack_out/Results.txt \
| wc -l
echo ""

echo "----------------------------------------------------------"

echo ""
echo "Number of loci _not_ matching miRBase miRNAs:"
awk '$21=="NA" {print $0}' ${output_dir_top}/ShortStack_out/Results.txt \
| wc -l

Number of loci matching miRBase miRNAs:
36

----------------------------------------------------------

Number of loci _not_ matching miRBase miRNAs:
19280

Although there are loci with matches to miRBase miRNAs, ShortStack did not annotate these clusters as miRNAs likely because they do not also match secondary structure criteria.

6.2.1 Directory tree of all ShortStack outputs

Many of these are large (by GitHub standards) BAM files, so will not be added to the repo.

Additionally, it’s unlikely we’ll utilize most of the other files (bigwig) generated by ShortStack.

# Load bash variables into memory
source .bashvars

tree -h ${output_dir_top}/

/home/shared/8TB_HDD_01/sam/gitrepos/deep-dive/D-Apul/output/13.1-Apul-sRNAseq-ShortStack-R1-reads/
├── [ 21K]  shortstack.log
└── [ 36K]  ShortStack_out
    ├── [ 28K]  alignment_details.tsv
    ├── [1.1M]  Counts.txt
    ├── [ 87K]  known_miRNAs.gff3
    ├── [1.8M]  known_miRNAs_unaligned.fasta
    ├── [5.3M]  merged_alignments_21_m.bw
    ├── [5.7M]  merged_alignments_21_p.bw
    ├── [5.0M]  merged_alignments_22_m.bw
    ├── [5.4M]  merged_alignments_22_p.bw
    ├── [ 10M]  merged_alignments_23-24_m.bw
    ├── [ 11M]  merged_alignments_23-24_p.bw
    ├── [1.4G]  merged_alignments.bam
    ├── [225K]  merged_alignments.bam.csi
    ├── [ 72M]  merged_alignments_other_m.bw
    ├── [ 75M]  merged_alignments_other_p.bw
    ├── [ 48M]  merged_alignments_sRNA-ACR-140-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bw
    ├── [ 54M]  merged_alignments_sRNA-ACR-145-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bw
    ├── [ 51M]  merged_alignments_sRNA-ACR-150-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bw
    ├── [ 44M]  merged_alignments_sRNA-ACR-173-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bw
    ├── [ 45M]  merged_alignments_sRNA-ACR-178-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bw
    ├── [7.8K]  mir.fasta
    ├── [1.9M]  Results.gff3
    ├── [2.9M]  Results.txt
    ├── [261M]  sRNA-ACR-140-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam
    ├── [227K]  sRNA-ACR-140-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam.csi
    ├── [299M]  sRNA-ACR-145-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam
    ├── [228K]  sRNA-ACR-145-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam.csi
    ├── [311M]  sRNA-ACR-150-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam
    ├── [231K]  sRNA-ACR-150-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam.csi
    ├── [275M]  sRNA-ACR-173-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam
    ├── [229K]  sRNA-ACR-173-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam.csi
    ├── [248M]  sRNA-ACR-178-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam
    ├── [230K]  sRNA-ACR-178-S1-TP2_R1_001.fastp-R1-31bp-auto_adapters-polyG.bam.csi
    └── [4.0K]  strucVis
        ├── [ 12K]  Cluster_10345.ps
        ├── [ 34K]  Cluster_10345.txt
        ├── [ 11K]  Cluster_10486.ps
        ├── [5.8K]  Cluster_10486.txt
        ├── [ 12K]  Cluster_10788.ps
        ├── [ 11K]  Cluster_10788.txt
        ├── [ 11K]  Cluster_10983.ps
        ├── [7.5K]  Cluster_10983.txt
        ├── [ 12K]  Cluster_1936.ps
        ├── [9.8K]  Cluster_1936.txt
        ├── [ 12K]  Cluster_2047.ps
        ├── [ 36K]  Cluster_2047.txt
        ├── [ 12K]  Cluster_2124.ps
        ├── [ 45K]  Cluster_2124.txt
        ├── [ 12K]  Cluster_2599.ps
        ├── [ 45K]  Cluster_2599.txt
        ├── [ 11K]  Cluster_2600.ps
        ├── [ 31K]  Cluster_2600.txt
        ├── [ 11K]  Cluster_3125.ps
        ├── [5.1K]  Cluster_3125.txt
        ├── [ 12K]  Cluster_3197.ps
        ├── [7.4K]  Cluster_3197.txt
        ├── [ 11K]  Cluster_323.ps
        ├── [ 20K]  Cluster_323.txt
        ├── [ 12K]  Cluster_3805.ps
        ├── [ 17K]  Cluster_3805.txt
        ├── [ 12K]  Cluster_4407.ps
        ├── [ 31K]  Cluster_4407.txt
        ├── [ 11K]  Cluster_5031.ps
        ├── [ 16K]  Cluster_5031.txt
        ├── [ 11K]  Cluster_532.ps
        ├── [3.7K]  Cluster_532.txt
        ├── [ 12K]  Cluster_5424.ps
        ├── [ 65K]  Cluster_5424.txt
        ├── [ 12K]  Cluster_568.ps
        ├── [ 55K]  Cluster_568.txt
        ├── [ 12K]  Cluster_6570.ps
        ├── [ 22K]  Cluster_6570.txt
        ├── [ 11K]  Cluster_6643.ps
        ├── [ 46K]  Cluster_6643.txt
        ├── [ 12K]  Cluster_7171.ps
        ├── [ 47K]  Cluster_7171.txt
        ├── [ 12K]  Cluster_7178.ps
        ├── [ 13K]  Cluster_7178.txt
        ├── [ 12K]  Cluster_7190.ps
        ├── [ 80K]  Cluster_7190.txt
        ├── [ 11K]  Cluster_7310.ps
        ├── [ 26K]  Cluster_7310.txt
        ├── [ 11K]  Cluster_8733.ps
        ├── [4.7K]  Cluster_8733.txt
        ├── [ 12K]  Cluster_8805.ps
        ├── [ 30K]  Cluster_8805.txt
        ├── [ 12K]  Cluster_9742.ps
        └── [ 31K]  Cluster_9742.txt

2 directories, 87 files

---

Citations

Axtell, Michael J. 2013. “ShortStack: Comprehensive Annotation and Quantification of Small RNA Genes.” RNA 19 (6): 740–51. https://doi.org/10.1261/rna.035279.112.

Johnson, Nathan R, Jonathan M Yeoh, Ceyda Coruh, and Michael J Axtell. 2016. “Improved Placement of Multi-Mapping Small RNAs.” G3 Genes|Genomes|Genetics 6 (7): 2103–11. https://doi.org/10.1534/g3.116.030452.

Shahid, Saima, and Michael J. Axtell. 2014. “Identification and Annotation of Small RNA Genes Using ShortStack.” Methods 67 (1): 20–27. https://doi.org/10.1016/j.ymeth.2013.10.004.