###
# Kraken

# Run Kraken
kraken2     --db=/scratch/sam/kraken2-db/     --paired     --threads 16     --report Calf_44.kraken2.report     Calf_44.R1.fastq Calf_44.R2.fastq > Calf_44.kraken2

# Run Bracken
bracken     -d /pub42/sam2/software/pipelines/db/bracken/     -i Calf_44.kraken2.report     -o Calf_44.bracken     -r 100     -l S     -t 8
mv Calf_44.kraken2_bracken.report Calf_44.bracken.tsv

# Filter
clamp_bracken_to_taxa.py     --input Calf_44.bracken.tsv     --rank S     --host_tax 9913,9606     --min_percent 0.1     --sample_id Calf_44     --output Calf_44.processed.tsv

###
# Humann (and MetaPhlAn)

# Run Humann
humann3      --input Calf_44.fastq     --output .     --search-mode uniref50     --o-log ./log     --threads 8     --metaphlan-options "-x v30_CHOCOPhlAn_201901 --bowtie2db /pub42/sam2/software/pipelines/db/metaphlan/choc/ -t rel_ab"     --diamond-options "--block-size 6.0 -c 1"
cp Calf_44_genefamilies.tsv Calf_44.genefamilies.tsv
cp Calf_44_pathabundance.tsv Calf_44.pathabundance.tsv
cp Calf_44_pathcoverage.tsv Calf_44.pathcoverage.tsv
cp Calf_44_humann_temp/Calf_44_metaphlan_bugs_list.tsv Calf_44.metaphlan.tsv

# Trim metaphlan columns
{ echo -e "clade\tCalf_44"; cat Calf_44.metaphlan.tsv | grep -v '^#' | cut -f1,3 ; } > Calf_44.metaphlan.trim.tsv

# Merge all metaphlan table files into one. Get all *.metaphlan.tsv files in the same directory and combine them into a single table
humann_join_tables     -i .     --file_name metaphlan     -o metaphlan.tsv

# Clamp metahplan table to only include species classifications
cat metaphlan.tsv | egrep '^clade|s__' | grep -v t__  | sed 's/.*|s__//g' > metaphlan.species.tsv

# Renormalise humann pathabundance and gene_family tables to cpm
humann_renorm_table     -u cpm     -i Calf_44.pathabundance.tsv     -s n     -o Calf_44.pathabundance.cpm.tsv
humann_renorm_table     -u cpm     -i Calf_44.genefamilies.tsv     -s n     -o Calf_44.genefamilies.cpm.tsv

# Convert gene_family tables to GOslim terms
bn=$(echo "Calf_44.genefamilies.cpm.tsv" | cut -f1 -d '.')
group_humann2_uniref_abundances_to_GO.sh     -i Calf_44.genefamilies.cpm.tsv     -m Calf_44.genefamilies.molecular_functions.tsv     -b Calf_44.genefamilies.biological_processes.tsv     -c Calf_44.genefamilies.cellular_components.tsv     -u /pub42/sam2/software/pipelines/db/go/map_go_uniref50.txt     -o ./

# Re-CPM
cut -f2,3 Calf_44.cellular_components.tsv > Calf_44.cellular_components.trim.tsv
humann_renorm_table     -u cpm     -i Calf_44.cellular_components.trim.tsv     -o Calf_44.cellular_components.cpm.tsv