GCAtools_prepare1.sh

#!/usr/bin/env bash

#     GCAtools
#     Copyright (C) Yumin Huang
#     Contact: ymhuang@fafu.edu.cn
# 
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.

gettime() {
        echo -e `date '+%Y-%m-%d %H:%M:%S ... '`
}

set -e

if [ -z "$1" ]; then
        echo "The arguments is empty!"
        exit 1
fi
#
Usage () {
	echo ""
	echo "Tool:  GCAtools prepare"
        echo "Usage: GCAtools_prepare -g <gff3 file> -f <genome file> [options]"
	echo ""
        echo "Description:"
        echo "  Prepare gene clusters for input"
        echo ""
	echo "Options:"
        echo "  -h Show this message and exit"
        echo "  -g <string> gff3 file"
        echo "  -f <string> genome fasta file"
        echo "  -c <string> cds fasta file [default False]"
        echo "  -p <string> pep fasta file [default False]"
        echo "  -d <string> output dir [default ./results]"
        echo "  -o <string> output prefix [default out]"
        echo ""
        echo "  -m <int> minimal cluster size [default 2]"
	echo ""
	echo "Author: Yumin Huang; ymhuang@fafu.edu.cn"
	echo ""
        exit 1
}
while getopts "hg:f:c:p:d:o:m:" opt
do
    case $opt in
        h)
                Usage
                exit 1
                ;;
        g)
                aname=$OPTARG
                ;;
        f)
                bname=$OPTARG
                ;;
        c)
                cname=$OPTARG
                ;;
        p)
                dname=$OPTARG
                ;;
        d)
                ename=$OPTARG
                ;;
        o)
                oname=$OPTARG
                ;;
        m)
                cls=$OPTARG
                ;;
        ?)
                echo "Unknow argument!"
                exit 1
                ;;
        esac
done

logo () {
        echo ""
        echo "   ==============================="
        echo "  ||                             ||"
        echo "  ||         GCAtools            ||"
        echo "  ||       Version: v0.0.1       ||"
        echo "  ||                             ||"
        echo "   ==============================="
        echo ""
}
logo

script_dir=$(dirname "$(realpath "${BASH_SOURCE[0]}")")
based=$(echo ${ename-./results})
size=$(echo ${cls-2})

if [[ "$size" =~ ^[0-9]+$ ]] && [ "$size" -gt 0 ] 2>/dev/null; then
    echo "Set minimal cluster size to ${size}..."
else
    echo "$size is not a positive integer..."
    exit;
fi

cds=$(echo ${cname-false})
pep=$(echo ${dname-false})

echo `gettime`"prepare files..."

rm -rf $based
mkdir -p $based

if [ "$cds" = false ]; then
    echo "no cds file input, generate cds file by gffread..."
    gffread $aname -g $bname -x ${based}/${oname}.cds.fa
    sed -i 's/>transcript:/>/g' ${based}/${oname}.cds.fa
    cds=$(echo ${based}/${oname}.cds.fa)
fi

if [ "$pep" = false ]; then
    echo "no pep file input, generate pep file by gffread..."
    gffread $aname -g $bname -y ${based}/${oname}.pep.fa
    sed -i 's/>transcript:/>/g' ${based}/${oname}.pep.fa
    pep=$(echo ${based}/${oname}.pep.fa)
fi

echo `gettime`"generate primary only dna/cds/pep files..."

python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id --primary_only $aname > ${based}/${oname}.bed
bioawk -c fastx '{print $1"\t"$2}' $pep | fgrep -w -f <(cut -f4 ${based}/${oname}.bed) - | awk '{print ">"$1"\n"$2}' > ${based}/${oname}.pep
bioawk -c fastx '{print $1"\t"$2}' $cds | fgrep -w -f <(cut -f4 ${based}/${oname}.bed) - | awk '{print ">"$1"\n"$2}' > ${based}/${oname}.cds
bedtools getfasta -fi $bname -bed ${based}/${oname}.bed -nameOnly -s | sed 's/(+)//g' | sed 's/(-)//g' > ${based}/${oname}.dna

echo `gettime`"run mmseqs..."

mmseqs easy-cluster ${based}/${oname}.pep ${based}/${oname}.pep_cluster tmp --min-seq-id 0.4 > ${based}/${oname}.mmseq.logs

cut -f1 ${based}/${oname}.pep_cluster_cluster.tsv | sort | uniq -c | awk -v var="$size" '$1>=var{print $2}' | fgrep -w -f - ${based}/${oname}.pep_cluster_cluster.tsv > ${based}/${oname}.pep_cluster_cluster.nosingleton.tsv
cut -f1 ${based}/${oname}.pep_cluster_cluster.nosingleton.tsv | uniq > ${based}/${oname}.cluster_names.txt

echo `gettime`"extract $(cat "${based}"/"${oname}".cluster_names.txt |wc -l) clusters..."

mkdir ${based}/${oname}.clusters_cds
mkdir ${based}/${oname}.clusters_dna
mkdir ${based}/${oname}.clusters_pep


python3 ${script_dir}/scripts/extract_clusters.py ${based}/${oname}.cluster_names.txt ${based}/${oname}.pep_cluster_all_seqs.fasta ${based}/${oname}.clusters_pep/

python3 ${script_dir}/scripts/dna_cluster.py ${based}/${oname}.pep_cluster_cluster.tsv ${based}/${oname}.cds > ${based}/${oname}.cds_cluster_all_seqs.fasta
python3 ${script_dir}/scripts/dna_cluster.py ${based}/${oname}.pep_cluster_cluster.tsv ${based}/${oname}.dna > ${based}/${oname}.dna_cluster_all_seqs.fasta

python3 ${script_dir}/scripts/extract_clusters.py ${based}/${oname}.cluster_names.txt ${based}/${oname}.cds_cluster_all_seqs.fasta ${based}/${oname}.clusters_cds/
python3 ${script_dir}/scripts/extract_clusters.py ${based}/${oname}.cluster_names.txt ${based}/${oname}.dna_cluster_all_seqs.fasta ${based}/${oname}.clusters_dna/

echo `gettime`"done..."