# Some stats
seqkit stats 1*S_*trimmed_ss*.fq.gz
file                                     format  type  num_seqs    sum_len  min_len  avg_len  max_len
16S_ERR1018543_1_trimmed_ssualign.fq.gz  FASTQ   DNA     26,496  7,348,681      255    277.4      296
16S_ERR1018543_2_trimmed_ssualign.fq.gz  FASTQ   DNA     26,496  7,851,808      250    296.3      305
16S_ERR1018546_1_trimmed_ssualign.fq.gz  FASTQ   DNA     31,554  8,748,253      250    277.2      296
16S_ERR1018546_2_trimmed_ssualign.fq.gz  FASTQ   DNA     31,554  9,373,913      263    297.1      305
18S_ERR1018543_1_trimmed_ssualign.fq.gz  FASTQ   DNA     10,164  2,818,891      276    277.3      296
18S_ERR1018543_2_trimmed_ssualign.fq.gz  FASTQ   DNA     10,164  3,013,853      285    296.5      305
18S_ERR1018546_1_trimmed_ssualign.fq.gz  FASTQ   DNA      4,248  1,177,629      276    277.2      296
18S_ERR1018546_2_trimmed_ssualign.fq.gz  FASTQ   DNA      4,248  1,262,690      263    297.2      305

seqkit stats 1*S_*trimmed.fq.gz
file                            format  type  num_seqs    sum_len  min_len  avg_len  max_len
16S_ERR1018543_1_trimmed.fq.gz  FASTQ   DNA     30,135  8,357,879      255    277.3      296
16S_ERR1018543_2_trimmed.fq.gz  FASTQ   DNA     30,135  8,931,640      250    296.4      305
16S_ERR1018546_1_trimmed.fq.gz  FASTQ   DNA     33,381  9,254,697      250    277.2      296
16S_ERR1018546_2_trimmed.fq.gz  FASTQ   DNA     33,381  9,916,913      263    297.1      305
18S_ERR1018543_1_trimmed.fq.gz  FASTQ   DNA      6,525  1,809,693      276    277.3      294
18S_ERR1018543_2_trimmed.fq.gz  FASTQ   DNA      6,525  1,934,021      285    296.4      305
18S_ERR1018546_1_trimmed.fq.gz  FASTQ   DNA      2,421    671,185      276    277.2      291
18S_ERR1018546_2_trimmed.fq.gz  FASTQ   DNA      2,421    719,690      285    297.3      305

# Reformat fastq to fasta
reformat.sh in=ERR1018543_1_trimmed.fq.gz in2=ERR1018543_2_trimmed.fq.gz out=ERR1018543_1_trimmed.fasta out2=ERR1018543_2_trimmed.fasta
reformat.sh in=ERR1018546_1_trimmed.fq.gz in2=ERR1018546_2_trimmed.fq.gz out=ERR1018546_1_trimmed.fasta out2=ERR1018546_2_trimmed.fasta

# Prepare ssu-align to be run in parallel
ssu-prep -f -x ERR1018543_1_trimmed.fasta ssu-align_ERR1018543_1 32
ssu-prep -f -x ERR1018543_2_trimmed.fasta ssu-align_ERR1018543_2 32
ssu-prep -f -x ERR1018546_1_trimmed.fasta ssu-align_ERR1018546_1 32
ssu-prep -f -x ERR1018546_2_trimmed.fasta ssu-align_ERR1018546_2 32

# Run ssu-align
./ssu-align_ERR1018543_1.ssu-align.sh
./ssu-align_ERR1018543_2.ssu-align.sh
./ssu-align_ERR1018546_1.ssu-align.sh
./ssu-align_ERR1018546_2.ssu-align.sh

# Get read ids for the 16S and 18S reads identified by ssu-align
PAIR=ERR1018543_1; find ssu-align_"${PAIR}" -regex '.*\(archaea.fa\|bacteria.fa\)$' -exec grep '>' {} \; | tr -d '>' | cut -f 1 -d ' ' > "${PAIR}"_16S_ssu-align.ids
PAIR=ERR1018543_2; find ssu-align_"${PAIR}" -regex '.*\(archaea.fa\|bacteria.fa\)$' -exec grep '>' {} \; | tr -d '>' | cut -f 1 -d ' ' > "${PAIR}"_16S_ssu-align.ids

PAIR=ERR1018546_1; find ssu-align_"${PAIR}" -regex '.*\(archaea.fa\|bacteria.fa\)$' -exec grep '>' {} \; | tr -d '>' | cut -f 1 -d ' ' > "${PAIR}"_16S_ssu-align.ids
PAIR=ERR1018546_2; find ssu-align_"${PAIR}" -regex '.*\(archaea.fa\|bacteria.fa\)$' -exec grep '>' {} \; | tr -d '>' | cut -f 1 -d ' ' > "${PAIR}"_16S_ssu-align.ids

PAIR=ERR1018543_1; find ssu-align_"${PAIR}" -regex '.*\(eukarya.fa\)$' -exec grep '>' {} \; | tr -d '>' | cut -f 1 -d ' ' > "${PAIR}"_18S_ssu-align.ids
PAIR=ERR1018543_2; find ssu-align_"${PAIR}" -regex '.*\(eukarya.fa\)$' -exec grep '>' {} \; | tr -d '>' | cut -f 1 -d ' ' > "${PAIR}"_18S_ssu-align.ids

PAIR=ERR1018546_1; find ssu-align_"${PAIR}" -regex '.*\(eukarya.fa\)$' -exec grep '>' {} \; | tr -d '>' | cut -f 1 -d ' ' > "${PAIR}"_18S_ssu-align.ids
PAIR=ERR1018546_2; find ssu-align_"${PAIR}" -regex '.*\(eukarya.fa\)$' -exec grep '>' {} \; | tr -d '>' | cut -f 1 -d ' ' > "${PAIR}"_18S_ssu-align.ids

# Split reads (we will only use a combination of FWD/REV ids, a more strict approach would be to select only those reads that have been in aligned in FWD/REV)
SSU=16S; SAMPLE=ERR1018543; filterbyname.sh in="${SAMPLE}"_1_trimmed.fq.gz in2="${SAMPLE}"_2_trimmed.fq.gz out="${SSU}"_"${SAMPLE}"_1_trimmed_ssualign.fq.gz out2="${SSU}"_"${SAMPLE}"_2_trimmed_ssualign.fq.gz names=<(cat "${SAMPLE}"_1_18S_ssu-align.ids "${SAMPLE}"_2_18S_ssu-align.ids) include=f
SSU=18S; SAMPLE=ERR1018543; filterbyname.sh in="${SAMPLE}"_1_trimmed.fq.gz in2="${SAMPLE}"_2_trimmed.fq.gz out="${SSU}"_"${SAMPLE}"_1_trimmed_ssualign.fq.gz out2="${SSU}"_"${SAMPLE}"_2_trimmed_ssualign.fq.gz names=<(cat "${SAMPLE}"_1_18S_ssu-align.ids "${SAMPLE}"_2_18S_ssu-align.ids) include=t

SSU=16S; SAMPLE=ERR1018546; filterbyname.sh in="${SAMPLE}"_1_trimmed.fq.gz in2="${SAMPLE}"_2_trimmed.fq.gz out="${SSU}"_"${SAMPLE}"_1_trimmed_ssualign.fq.gz out2="${SSU}"_"${SAMPLE}"_2_trimmed_ssualign.fq.gz names=<(cat "${SAMPLE}"_1_"${SSU}"_ssu-align.ids "${SAMPLE}"_2_"${SSU}"_ssu-align.ids) include=f
SSU=18S; SAMPLE=ERR1018546; filterbyname.sh in="${SAMPLE}"_1_trimmed.fq.gz in2="${SAMPLE}"_2_trimmed.fq.gz out="${SSU}"_"${SAMPLE}"_1_trimmed_ssualign.fq.gz out2="${SSU}"_"${SAMPLE}"_2_trimmed_ssualign.fq.gz names=<(cat "${SAMPLE}"_1_18S_ssu-align.ids "${SAMPLE}"_2_18S_ssu-align.ids) include=t