Use miRDeep2 (Friedländer et al. 2011) to identify potential miRNAs using P.meandrina sRNAseq reads. The P.meandrina genome will be used as the reference genome. A modified version of the mature miRBase miRNA database, which includes cnidarian miRNA culled from literature by Jill Ahsley, will also be used.
Inputs:
Requires collapsed reads (i.e. concatenated, unique reads) in FastA format. See 10.1-Pmea-sRNAseq-BLASTn-31bp-fastp-merged-cnidarian_miRBase.Rmd for code.
P.meandrina genome FastA. See 12-Pmea-sRNAseq-MirMachine.Rmd for download info if needed.
Utilizes a modified version, which includes cnidarian miRNA culled from literature by Jill Ahsley.
Outputs:
Primary outputs are a result table in BED, CSV (tab-delimited), and HTML formats.
Due to the nature of mirDeep2’s naming, trying to use variable
names is challenging. As such, the chunks processing those files will
require manual intervention to identify and provide the output
filename(s); they are not handled at in the
.bashvars file at the top of this script.
Other output files are too large for GitHub (and some (all?) are not needed for this analysis). Please find a full backup here:
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-adapters-polyG-31bp-merged.fq.gz'"
echo ""
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}/F-Pmea/output/11.1-Pmea-sRNAseq-miRdeep2-31bp-fastp-merged-cnidarian_miRBase'
echo 'export genome_fasta_dir=${deep_dive_dir}/F-Pmea/data'
echo 'export trimmed_fastqs_dir="${deep_dive_dir}/F-Pmea/output//08.2-Pmea-sRNAseq-trimming-31bp-fastp-merged/trimmed-reads"'
echo 'export collapsed_reads_dir="${deep_dive_dir}/F-Pmea/output/10.1-Pmea-sRNAseq-BLASTn-31bp-fastp-merged-cnidarian_miRBase"'
echo ""
echo "# Input/Output files"
echo 'export collapsed_reads_fasta="collapsed-reads-all.fasta"'
echo 'export collapsed_reads_mirdeep2="collapsed-reads-all-mirdeep2.fasta"'
echo 'export collapsed_reads_gt17bp_mirdeep2="collapsed_reads_gt17bp_mirdeep2.fasta"'
echo 'export concatenated_trimmed_reads_fastq="concatenated-trimmed-reads-all.fastq.gz"'
echo 'export genome_fasta_name="Pocillopora_meandrina_HIv1.assembly.fasta"'
echo 'export genome_fasta_no_spaces="Pocillopora_meandrina_HIv1_nospaces.assembly.fasta"'
echo 'export mirdeep2_mapping_file="Pmea-mirdeep2-mapping.arf"'
echo 'export mirbase_mature_fasta_name="cnidarian-mirbase-mature-v22.1.fasta"'
echo 'export mirbase_mature_fasta_no_spaces="cnidarian-mirbase-mature-v22.1-no_spaces.fa"'
echo ""
echo "# Paths to programs"
echo 'export mirdeep2_mapper="mapper.pl"'
echo 'export mirdeep2="miRDeep2.pl"'
echo 'export mirdeep2_fastaparse="fastaparse.pl"'
echo 'export bowtie_build="/home/shared/bowtie-1.3.1-linux-x86_64/bowtie-build"'
echo ""
echo "# Set number of CPUs to use"
echo 'export threads=46'
echo ""
echo "# Initialize arrays"
echo 'export trimmed_fastqs_array=()'
} > .bashvars
cat .bashvars#### Assign Variables ####
# Trimmed FastQ naming pattern
export trimmed_fastqs_pattern='*fastp-adapters-polyG-31bp-merged.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}/F-Pmea/output/11.1-Pmea-sRNAseq-miRdeep2-31bp-fastp-merged-cnidarian_miRBase
export genome_fasta_dir=${deep_dive_dir}/F-Pmea/data
export trimmed_fastqs_dir="${deep_dive_dir}/F-Pmea/output//08.2-Pmea-sRNAseq-trimming-31bp-fastp-merged/trimmed-reads"
export collapsed_reads_dir="${deep_dive_dir}/F-Pmea/output/10.1-Pmea-sRNAseq-BLASTn-31bp-fastp-merged-cnidarian_miRBase"
# Input/Output files
export collapsed_reads_fasta="collapsed-reads-all.fasta"
export collapsed_reads_mirdeep2="collapsed-reads-all-mirdeep2.fasta"
export collapsed_reads_gt17bp_mirdeep2="collapsed_reads_gt17bp_mirdeep2.fasta"
export concatenated_trimmed_reads_fastq="concatenated-trimmed-reads-all.fastq.gz"
export genome_fasta_name="Pocillopora_meandrina_HIv1.assembly.fasta"
export genome_fasta_no_spaces="Pocillopora_meandrina_HIv1_nospaces.assembly.fasta"
export mirdeep2_mapping_file="Pmea-mirdeep2-mapping.arf"
export mirbase_mature_fasta_name="cnidarian-mirbase-mature-v22.1.fasta"
export mirbase_mature_fasta_no_spaces="cnidarian-mirbase-mature-v22.1-no_spaces.fa"
# Paths to programs
export mirdeep2_mapper="mapper.pl"
export mirdeep2="miRDeep2.pl"
export mirdeep2_fastaparse="fastaparse.pl"
export bowtie_build="/home/shared/bowtie-1.3.1-linux-x86_64/bowtie-build"
# Set number of CPUs to use
export threads=46
# Initialize arrays
export trimmed_fastqs_array=()Per miRDeep2 documentation:
The readID must end with _xNumber and is not allowed to contain whitespaces. has to have the format name_uniqueNumber_xnumber
# Load bash variables into memory
source .bashvars
# Append miRDeep2 to system PATH and set PERL5LIB
export PATH=$PATH:/home/shared/mirdeep2/bin
export PERL5LIB=$PERL5LIB:/home/shared/mirdeep2/lib/perl5
mkdir --parents "${output_dir_top}"
sed '/^>/ s/-/_x/g' "${collapsed_reads_dir}/${collapsed_reads_fasta}" \
| sed '/^>/ s/>/>seq_/' \
> "${output_dir_top}/${collapsed_reads_mirdeep2}"
# Filter for reads at least 17 bases long
# Min. read length required for MirDeep2
${mirdeep2_fastaparse} \
"${output_dir_top}/${collapsed_reads_mirdeep2}" \
-a 17 \
> "${output_dir_top}/${collapsed_reads_gt17bp_mirdeep2}" \
2> "${output_dir_top}/fasta_parse.log"
grep "^>" ${collapsed_reads_dir}/${collapsed_reads_fasta} \
| head
echo ""
echo "--------------------------------------------------"
echo ""
grep "^>" "${output_dir_top}/${collapsed_reads_mirdeep2}" \
| head
echo ""
echo "--------------------------------------------------"
echo ""
grep "^>" ${output_dir_top}/${collapsed_reads_gt17bp_mirdeep2} \
| head>1-10257663
>2-2952210
>3-1527059
>4-1205244
>5-1189262
>6-630167
>7-548099
>8-512116
>9-432671
>10-407718
--------------------------------------------------
>seq_1_x10257663
>seq_2_x2952210
>seq_3_x1527059
>seq_4_x1205244
>seq_5_x1189262
>seq_6_x630167
>seq_7_x548099
>seq_8_x512116
>seq_9_x432671
>seq_10_x407718
--------------------------------------------------
>seq_1_x10257663
>seq_2_x2952210
>seq_3_x1527059
>seq_4_x1205244
>seq_5_x1189262
>seq_6_x630167
>seq_7_x548099
>seq_8_x512116
>seq_9_x432671
>seq_10_x407718miRDeep2 can’t process genome FastAs with spaces in the description lines.
So, I’m replacing spaces with underscores.
And, for aesthetics, I’m also removing commas.
# Load bash variables into memory
source .bashvars
sed '/^>/ s/ /_/g' "${genome_fasta_dir}/${genome_fasta_name}" \
| sed '/^>/ s/,//g' \
> "${genome_fasta_dir}/${genome_fasta_no_spaces}"
grep "^>" ${genome_fasta_dir}/${genome_fasta_name} \
| head
echo ""
echo "--------------------------------------------------"
echo ""
grep "^>" ${genome_fasta_dir}/${genome_fasta_no_spaces} \
| head>Pocillopora_meandrina_HIv1___Sc0000000
>Pocillopora_meandrina_HIv1___Sc0000001
>Pocillopora_meandrina_HIv1___Sc0000002
>Pocillopora_meandrina_HIv1___Sc0000003
>Pocillopora_meandrina_HIv1___Sc0000004
>Pocillopora_meandrina_HIv1___Sc0000005
>Pocillopora_meandrina_HIv1___Sc0000006
>Pocillopora_meandrina_HIv1___Sc0000007
>Pocillopora_meandrina_HIv1___Sc0000008
>Pocillopora_meandrina_HIv1___Sc0000009
--------------------------------------------------
>Pocillopora_meandrina_HIv1___Sc0000000
>Pocillopora_meandrina_HIv1___Sc0000001
>Pocillopora_meandrina_HIv1___Sc0000002
>Pocillopora_meandrina_HIv1___Sc0000003
>Pocillopora_meandrina_HIv1___Sc0000004
>Pocillopora_meandrina_HIv1___Sc0000005
>Pocillopora_meandrina_HIv1___Sc0000006
>Pocillopora_meandrina_HIv1___Sc0000007
>Pocillopora_meandrina_HIv1___Sc0000008
>Pocillopora_meandrina_HIv1___Sc0000009miRDeep2 can’t process genome FastAs with spaces in the description lines.
So, I’m replacing spaces with underscores.
And, for aesthetics, I’m also removing commas.
# Load bash variables into memory
source .bashvars
sed '/^>/ s/ /_/g' "${deep_dive_data_dir}/${mirbase_mature_fasta_name}" \
| sed '/^>/ s/,//g' \
> "${deep_dive_data_dir}/${mirbase_mature_fasta_no_spaces}"
grep "^>" ${deep_dive_data_dir}/${mirbase_mature_fasta_name} \
| head
echo ""
echo "--------------------------------------------------"
echo ""
grep "^>" ${deep_dive_data_dir}/${mirbase_mature_fasta_no_spaces} \
| head>cel-let-7-5p MIMAT0000001 Caenorhabditis elegans let-7-5p
>cel-let-7-3p MIMAT0015091 Caenorhabditis elegans let-7-3p
>cel-lin-4-5p MIMAT0000002 Caenorhabditis elegans lin-4-5p
>cel-lin-4-3p MIMAT0015092 Caenorhabditis elegans lin-4-3p
>cel-miR-1-5p MIMAT0020301 Caenorhabditis elegans miR-1-5p
>cel-miR-1-3p MIMAT0000003 Caenorhabditis elegans miR-1-3p
>cel-miR-2-5p MIMAT0020302 Caenorhabditis elegans miR-2-5p
>cel-miR-2-3p MIMAT0000004 Caenorhabditis elegans miR-2-3p
>cel-miR-34-5p MIMAT0000005 Caenorhabditis elegans miR-34-5p
>cel-miR-34-3p MIMAT0015093 Caenorhabditis elegans miR-34-3p
--------------------------------------------------
>cel-let-7-5p_MIMAT0000001_Caenorhabditis_elegans_let-7-5p
>cel-let-7-3p_MIMAT0015091_Caenorhabditis_elegans_let-7-3p
>cel-lin-4-5p_MIMAT0000002_Caenorhabditis_elegans_lin-4-5p
>cel-lin-4-3p_MIMAT0015092_Caenorhabditis_elegans_lin-4-3p
>cel-miR-1-5p_MIMAT0020301_Caenorhabditis_elegans_miR-1-5p
>cel-miR-1-3p_MIMAT0000003_Caenorhabditis_elegans_miR-1-3p
>cel-miR-2-5p_MIMAT0020302_Caenorhabditis_elegans_miR-2-5p
>cel-miR-2-3p_MIMAT0000004_Caenorhabditis_elegans_miR-2-3p
>cel-miR-34-5p_MIMAT0000005_Caenorhabditis_elegans_miR-34-5p
>cel-miR-34-3p_MIMAT0015093_Caenorhabditis_elegans_miR-34-3pmiRDeep2 requires a Bowtie v1 genome index - cannot use Bowtie2 index
# Load bash variables into memory
source .bashvars
# Check for existence of genome index first
if [ ! -f "${genome_fasta_no_spaces%.*}.*ebwt" ]; then
${bowtie_build} \
${genome_fasta_dir}/${genome_fasta_no_spaces} \
${genome_fasta_dir}/${genome_fasta_no_spaces%.*} \
--threads ${threads} \
--quiet
fiRequires genome to be previously indexed with Bowtie.
Additionally, requires user to enter path to their
mirdeep2 directory as well as their perl5
installation.
# Load bash variables into memory
source .bashvars
# Append miRDeep2 to system PATH and set PERL5LIB
export PATH=$PATH:/home/shared/mirdeep2/bin
export PERL5LIB=$PERL5LIB:/home/shared/mirdeep2/lib/perl5
# Run miRDeep2 mapping
time \
${mirdeep2_mapper} \
${output_dir_top}/${collapsed_reads_gt17bp_mirdeep2} \
-c \
-p ${genome_fasta_dir}/${genome_fasta_no_spaces%.*} \
-t ${output_dir_top}/${mirdeep2_mapping_file} \
-o ${threads}Recommendation is to use the closest related species in the miRDeep2 options, even if the species isn’t very closely related. The documentation indicates that miRDeep2 is always more accurate when at least a species is provided.
The options provided to the command are as follows:
none: Known miRNAs of the species being analyzed.-t S.pupuratus: Related species.none: Known miRNA precursors in this species.-P: Specifies miRBase version > 18.-v: Remove temporary files after completion.-g -1: Number of precursors to anlayze. A setting of
-1 will analyze all. Default is 50,000 I set this to
-1 after multiple attempts to run using the default kept
failing.NOTE: This will take an extremely long time to run (days).
Could possible by shortened by excluding randfold
analysis.
# Load bash variables into memory
source .bashvars
# Append miRDeep2 to system PATH and set PERL5LIB
export PATH=$PATH:/home/shared/mirdeep2/bin
export PERL5LIB=$PERL5LIB:/home/shared/mirdeep2/lib/perl5
time \
${mirdeep2} \
${output_dir_top}/${collapsed_reads_gt17bp_mirdeep2} \
${genome_fasta_dir}/${genome_fasta_no_spaces} \
${output_dir_top}/${mirdeep2_mapping_file} \
none \
${deep_dive_data_dir}/${mirbase_mature_fasta_no_spaces} \
none \
-t S.purpuratus \
-P \
-v \
-g -1 \
2>${output_dir_top}/miRDeep2-S.purpuratus-report.log# Load bash variables into memory
source .bashvars
tail -n 6 ${output_dir_top}/miRDeep2-S.purpuratus-report.logmiRDeep runtime:
started: 15:27:22
ended: 4:31:50
total:37h:4m:28sMiRDeep2 outputs all files to the current working directly with no way to redirect so want to move to intended output directory.
Output files will be in the format of result_* and
error_*
# Load bash variables into memory
source .bashvars
for file in result_* error_*
do
mv "${file}" "${output_dir_top}/"
done# Load bash variables into memory
source .bashvars
ls -l | grep "^d"drwxr-xr-x 4 sam sam 4096 Feb 15 15:16 08.1-Pmea-sRNAseq-trimming-R1-only_cache
drwxr-xr-x 4 sam sam 4096 Feb 16 10:13 08.2-Pmea-sRNAseq-trimming-31bp-fastp-merged_cache
drwxr-xr-x 4 sam sam 4096 Apr 22 09:52 10.1-Pmea-sRNAseq-BLASTn-31bp-fastp-merged-cnidarian_miRBase_cache
drwxr-xr-x 4 sam sam 4096 Apr 19 12:24 10.2-Pmea-sRNAseq-BLASTn-31bp-fastp-merged_cache
drwxr-xr-x 4 sam sam 4096 Nov 15 13:47 10-Pmea-sRNAseq-BLASTn_cache
drwxr-xr-x 4 sam sam 4096 Feb 16 14:31 13.1.1-Pmea-sRNAseq-ShortStack-R1-reads-cnidarian_miRBase_cache
drwxr-xr-x 4 sam sam 4096 Feb 16 11:36 13.1-Pmea-sRNAseq-ShortStack-R1-reads_cache
drwxr-xr-x 4 sam sam 4096 Feb 16 16:51 13.2.1-Pmea-sRNAseq-ShortStack-31bp-fastp-merged-cnidarian_miRBase_cache
drwxr-xr-x 4 sam sam 4096 Feb 16 12:36 13.2-Pmea-sRNAseq-ShortStack-31bp-fastp-merged_cache
drwxr-xr-x 4 sam sam 4096 Dec 1 09:42 13-Pmea-sRNAseq-ShortStack_cache
drwxr-xr-x 4 sam sam 4096 Dec 1 09:42 analyses
drwxr-xr-x 3 sam sam 4096 Dec 1 09:42 data# Load bash variables into memory
source .bashvars
rsync -aP . --include='dir***' --exclude='*' --quiet "${output_dir_top}/"
rsync -aP . --include='map***' --exclude='*' --quiet "${output_dir_top}/"
rsync -aP . --exclude='mirgene*' --include='mir***' --exclude='*' --quiet "${output_dir_top}/"
rsync -aP . --include='pdfs***' --exclude='*' --quiet "${output_dir_top}/"
echo ""
echo "Check new location:"
echo ""
ls -l "${output_dir_top}/" | grep "^d"
FYI - eval=FALSE is set because the following command
will only work once…
ls --directory dir_* mapper_logs mirdeep_runs mirna_results* pdfs_*# Load bash variables into memory
source .bashvars
rm -rf dir_* mapper_logs mirdeep_runs mirna_results* pdfs_*
ls --directory dir_* mapper_logs mirdeep_runs mirna_results* pdfs_*ls: cannot access 'dir_*': No such file or directory
ls: cannot access 'mapper_logs': No such file or directory
ls: cannot access 'mirdeep_runs': No such file or directory
ls: cannot access 'mirna_results*': No such file or directory
ls: cannot access 'pdfs_*': No such file or directoryNOTE: You will need to replace the name of the output CSV with the one from the current run.
The formatting of this CSV is terrible. It has a 3-column table on top of a 17-column table. This makes parsing a bit of a pain in its raw format.
# Load bash variables into memory
source .bashvars
head -n 30 "${output_dir_top}/result_22_04_2024_t_15_27_22.csv"miRDeep2 score estimated signal-to-noise excision gearing
10 6.9 1
9 7.1 1
8 7.1 1
7 7 1
6 7.4 1
5 9.5 1
4 7.5 1
3 4.5 1
2 2.7 1
1 2.1 1
0 2.3 1
-1 2.4 1
-2 1.5 1
-3 0.8 1
-4 0.6 1
-5 0.5 1
-6 0.6 1
-7 0.6 1
-8 0.7 1
-9 0.7 1
-10 0.8 1
novel miRNAs predicted by miRDeep2
provisional id miRDeep2 score estimated probability that the miRNA candidate is a true positive rfam alert total read count mature read count loop read count star read count significant randfold p-value miRBase miRNA example miRBase miRNA with the same seed UCSC browser NCBI blastn consensus mature sequence consensus star sequence consensus precursor sequence precursor coordinate
Pocillopora_meandrina_HIv1___Sc0000000_60202 266239.8 - 522212 520540 0 1672 yes - mmu-miR-710_MIMAT0003500_Mus_musculus_miR-710 - - ucaagucuaggcugguuaguuu cuaaaccagacuaggcuucagc cuaaaccagacuaggcuucagcauauuuauuuugucaagucuaggcugguuaguuu Pocillopora_meandrina_HIv1___Sc0000000:20372434..20372490:-
Pocillopora_meandrina_HIv1___Sc0000017_655645 67675.4 - 132734 82659 0 50075 yes - tca-miR-6008-5p_MIMAT0023557_Tribolium_castaneum_miR-6008-5p - - aacugcugagauucuauggauuu auccauagaacuucugcauuuga aacugcugagauucuauggauuuaauuuaaauccauagaacuucugcauuuga Pocillopora_meandrina_HIv1___Sc0000017:5050927..5050980:-
Pocillopora_meandrina_HIv1___Sc0000005_273871 65863.2 - 129179 127776 0 1403 yes - nve-miR-2023-3p_MIMAT0009756_Nematostella_vectensis_miR-2023-3p - - aaagaaguacaagugguaggg cugccacucguauuuucuuuca cugccacucguauuuucuuucacguuuaucgaugaaagaaguacaagugguaggg Pocillopora_meandrina_HIv1___Sc0000005:601591..601646:+This will match the line beginning with provisional id
and print to the end of the file (represented by the $p.
$ = end, p = print)
# Load bash variables into memory
source .bashvars
sed --quiet '/provisional id/,$p' "${output_dir_top}/result_22_04_2024_t_15_27_22.csv" \
> "${output_dir_top}/parsable-result_22_04_2024_t_15_27_22.csv"
head "${output_dir_top}/parsable-result_22_04_2024_t_15_27_22.csv"provisional id miRDeep2 score estimated probability that the miRNA candidate is a true positive rfam alert total read count mature read count loop read count star read count significant randfold p-value miRBase miRNA example miRBase miRNA with the same seed UCSC browser NCBI blastn consensus mature sequence consensus star sequence consensus precursor sequence precursor coordinate
Pocillopora_meandrina_HIv1___Sc0000000_60202 266239.8 - 522212 520540 0 1672 yes - mmu-miR-710_MIMAT0003500_Mus_musculus_miR-710 - - ucaagucuaggcugguuaguuu cuaaaccagacuaggcuucagc cuaaaccagacuaggcuucagcauauuuauuuugucaagucuaggcugguuaguuu Pocillopora_meandrina_HIv1___Sc0000000:20372434..20372490:-
Pocillopora_meandrina_HIv1___Sc0000017_655645 67675.4 - 132734 82659 0 50075 yes - tca-miR-6008-5p_MIMAT0023557_Tribolium_castaneum_miR-6008-5p - - aacugcugagauucuauggauuu auccauagaacuucugcauuuga aacugcugagauucuauggauuuaauuuaaauccauagaacuucugcauuuga Pocillopora_meandrina_HIv1___Sc0000017:5050927..5050980:-
Pocillopora_meandrina_HIv1___Sc0000005_273871 65863.2 - 129179 127776 0 1403 yes - nve-miR-2023-3p_MIMAT0009756_Nematostella_vectensis_miR-2023-3p - - aaagaaguacaagugguaggg cugccacucguauuuucuuuca cugccacucguauuuucuuucacguuuaucgaugaaagaaguacaagugguaggg Pocillopora_meandrina_HIv1___Sc0000005:601591..601646:+
Pocillopora_meandrina_HIv1___xfSc0000017_985353 65601.9 - 128675 126663 30 1982 no - dme-miR-6-1-5p_MIMAT0020787_Drosophila_melanogaster_miR-6-1-5p - - ugggaauaccccguguuguaggcuu agcccugucgggcgggguuaguacuug agcccugucgggcgggguuaguacuugggaugggugaccgccugggaauaccccguguuguaggcuu Pocillopora_meandrina_HIv1___xfSc0000017:39556..39623:-
Pocillopora_meandrina_HIv1___Sc0000009_442963 64606.4 - 126723 126663 17 43 no - dme-miR-6-1-5p_MIMAT0020787_Drosophila_melanogaster_miR-6-1-5p - - ugggaauaccccguguuguaggcuu gcccugucgggcgggguu gcccugucgggcgggguuuaguacuuggaugggugaccgccugggaauaccccguguuguaggcuu Pocillopora_meandrina_HIv1___Sc0000009:11978171..11978237:-
Pocillopora_meandrina_HIv1___Sc0000002_166394 54159.5 - 106223 105945 0 278 yes - tca-miR-3860-3p_MIMAT0018793_Tribolium_castaneum_miR-3860-3p - - uuguguaacucccuaaggaagg cucuuucgggugucacacaacg cucuuucgggugucacacaacgucgucaaggagcguuguguaacucccuaaggaagg Pocillopora_meandrina_HIv1___Sc0000002:3841965..3842022:-
Pocillopora_meandrina_HIv1___Sc0000016_620672 45914.8 - 90052 84515 0 5537 yes - rlcv-miR-rL1-29-5p_MIMAT0019194_Rhesus_lymphocryptovirus_miR-rL1-29-5p - - ucagucccaccaucucaccaau ggugagcuguuuggacuuaua ggugagcuguuuggacuuauauuauugguaucagucccaccaucucaccaau Pocillopora_meandrina_HIv1___Sc0000016:7550614..7550666:+
Pocillopora_meandrina_HIv1___Sc0000016_620687 45794.3 - 89817 84515 0 5302 yes - rlcv-miR-rL1-29-5p_MIMAT0019194_Rhesus_lymphocryptovirus_miR-rL1-29-5p - - ucagucccaccaucucaccaau ggugagcuguuugaacuuaua ggugagcuguuugaacuuauauuauugguaucagucccaccaucucaccaau Pocillopora_meandrina_HIv1___Sc0000016:7551980..7552032:+
Pocillopora_meandrina_HIv1___Sc0000016_620701 45794.3 - 89817 84515 0 5302 yes - rlcv-miR-rL1-29-5p_MIMAT0019194_Rhesus_lymphocryptovirus_miR-rL1-29-5p - - ucagucccaccaucucaccaau ggugagcuguuugaacuuaua ggugagcuguuugaacuuauauuauugguaucagucccaccaucucaccaau Pocillopora_meandrina_HIv1___Sc0000016:7553123..7553175:+This chunk provides a more concise overview of the data and it’s columns.
mirdeep_result.df <- read.csv("../output/11.1-Pmea-sRNAseq-miRdeep2-31bp-fastp-merged-cnidarian_miRBase/parsable-result_22_04_2024_t_15_27_22.csv",
header = TRUE,
sep = "\t")
str(mirdeep_result.df)'data.frame': 1040 obs. of 17 variables:
$ provisional.id : chr "Pocillopora_meandrina_HIv1___Sc0000000_60202" "Pocillopora_meandrina_HIv1___Sc0000017_655645" "Pocillopora_meandrina_HIv1___Sc0000005_273871" "Pocillopora_meandrina_HIv1___xfSc0000017_985353" ...
$ miRDeep2.score : num 266240 67675 65863 65602 64606 ...
$ estimated.probability.that.the.miRNA.candidate.is.a.true.positive: logi NA NA NA NA NA NA ...
$ rfam.alert : chr "-" "-" "-" "-" ...
$ total.read.count : int 522212 132734 129179 128675 126723 106223 90052 89817 89817 89808 ...
$ mature.read.count : int 520540 82659 127776 126663 126663 105945 84515 84515 84515 84515 ...
$ loop.read.count : int 0 0 0 30 17 0 0 0 0 0 ...
$ star.read.count : int 1672 50075 1403 1982 43 278 5537 5302 5302 5293 ...
$ significant.randfold.p.value : chr "yes" "yes" "yes" "no" ...
$ miRBase.miRNA : chr "-" "-" "-" "-" ...
$ example.miRBase.miRNA.with.the.same.seed : chr "mmu-miR-710_MIMAT0003500_Mus_musculus_miR-710" "tca-miR-6008-5p_MIMAT0023557_Tribolium_castaneum_miR-6008-5p" "nve-miR-2023-3p_MIMAT0009756_Nematostella_vectensis_miR-2023-3p" "dme-miR-6-1-5p_MIMAT0020787_Drosophila_melanogaster_miR-6-1-5p" ...
$ UCSC.browser : chr "-" "-" "-" "-" ...
$ NCBI.blastn : chr "-" "-" "-" "-" ...
$ consensus.mature.sequence : chr "ucaagucuaggcugguuaguuu" "aacugcugagauucuauggauuu" "aaagaaguacaagugguaggg" "ugggaauaccccguguuguaggcuu" ...
$ consensus.star.sequence : chr "cuaaaccagacuaggcuucagc" "auccauagaacuucugcauuuga" "cugccacucguauuuucuuuca" "agcccugucgggcgggguuaguacuug" ...
$ consensus.precursor.sequence : chr "cuaaaccagacuaggcuucagcauauuuauuuugucaagucuaggcugguuaguuu" "aacugcugagauucuauggauuuaauuuaaauccauagaacuucugcauuuga" "cugccacucguauuuucuuucacguuuaucgaugaaagaaguacaagugguaggg" "agcccugucgggcgggguuaguacuugggaugggugaccgccugggaauaccccguguuguaggcuu" ...
$ precursor.coordinate : chr "Pocillopora_meandrina_HIv1___Sc0000000:20372434..20372490:-" "Pocillopora_meandrina_HIv1___Sc0000017:5050927..5050980:-" "Pocillopora_meandrina_HIv1___Sc0000005:601591..601646:+" "Pocillopora_meandrina_HIv1___xfSc0000017:39556..39623:-" ...This provides some rudimentary numbers for the miRDeep2 output.
Further analysis is possibly desired to evaluate score thresholds, miRNA families, etc.
# Load bash variables into memory
source .bashvars
# Total predicted miRNAS
total_miRNAs=$(awk 'NR > 1' ${output_dir_top}/parsable-result_22_04_2024_t_15_27_22.csv \
| wc -l
)
echo "Total of predicted miRNAs: ${total_miRNAs}"
echo ""
# Matches to known mature miRNAs
mature_miRNAs=$(awk -F'\t' '$11 != "-" && $11 != "" {print $11}' ${output_dir_top}/parsable-result_22_04_2024_t_15_27_22.csv \
| wc -l
)
echo "Number of seed matches to known miRNAS: ${mature_miRNAs}"
echo ""
# Novel miRNAs
novel_miRNAs=$(awk -F "\t" '$11 == "-" || $11 == "" {print $11}' ${output_dir_top}/parsable-result_22_04_2024_t_15_27_22.csv \
| awk 'NR > 1' \
| wc -l
)
echo "Number of novel miRNAs: ${novel_miRNAs}"
Total of predicted miRNAs: 1040
Number of seed matches to known miRNAS: 912
Number of novel miRNAs: 128