#! /bin/bash
echo "Before you run this script make sure you unzip the ssu and refseq fasta files. Then we edit the GTDB ids to replace the . to a _"
cat ssu_all_r202.fna | sed -e 's/\.1/_1/g' -e 's/\.2/_2/g' -e 's/\.3/_3/g' > newSSU_all_r202.fa
echo "The following script will rename the GTDB fasta files to suit DADA2 format"

cat newSSU_all_r202.fa | awk -F " " '/^>/{print $1" "$2" "$3; next}{print $1}' > all_r202.fa

cat all_r202.fa | awk -F "~" '/^>/{print ">"$2" "$3"("$1")"; next}{print $1}' > GTDB_ssu_all_r202.fa

cat GTDB_ssu_all_r202.fa | awk -F " " '/^>/{print ">"$2" "$3""$4; next}{print $1}' > aGTDB_ssu_all_r202.fa


sed -e 's/d__//g' -e 's/p__//g' -e 's/c__//g' -e 's/o__//g' -e 's/f__//g' -e 's/g__//g' -e 's/s__//g' -e 's/(>/(/g' bGTDB_ssu_all_r202.fa > all_GTDBr202_16S_seqs.fa

echo "Number of all sequences"
grep -c '>' all_GTDBr202_16S_seqs.fa

echo "Number of bacterial sequences"
grep -c 'Bacteria' all_GTDBr202_16S_seqs.fa

echo "Number of archaeal sequences"
grep -c 'Archaea' all_GTDBr202_16S_seqs.fa

echo "Now we extract ATCC even genomic mock community species 16S gene using GTDB ids (the ids were extyracted from the GTDB_metadata file by searching genbank ids"

seqkit grep -r --by-name --pattern-file ATCC_GTDBr202_IDs.txt all_GTDBr202_16S_seqs.fa --line-width 0 > ATCC_Even_Mock_community_GTDBr202.fa

echo "you can limit the sequence lengths as follows"
#seqkit seq -m 500 -u -w 0 -j 4 -g ATCC_Even_Mock_community_GTDBr202.fa > reducedATCC_Even_Mock_community_GTDBr202.fa

echo "The follwing number is how many sequences match to ATCC 1002 in the GTDB release 202"
grep -c "^>" ATCC_Even_Mock_community_GTDBr202.fa

rm newSSU_all_r202.fa
rm all_r202.fa
rm GTDB_ssu_all_r202.fa
rm aGTDB_ssu_all_r202.fa
rm bGTDB_ssu_all_r202.fa

echo "Now we do the same for RefSeq but this time we use the zenodo uploaded database (less primitive code)"

seqkit grep -r --by-name --pattern-file mock_species.txt RefSeq_16S_6-11-20_RDPv16_fullTaxo.fa --line-width 0 > ATCC_Even_Mock_community_RefSeq.fa

grep -c "^>" ATCC_Even_Mock_community_RefSeq.fa

echo "THE END"
