-
Notifications
You must be signed in to change notification settings - Fork 5
/
FunOMIC.sh
executable file
·156 lines (119 loc) · 6.32 KB
/
FunOMIC.sh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
#!/usr/bin/env bash
################################################################################
Help()
{
# Display Help
echo
echo "FunOMIC: mycobiome taxonomic and functional profiling pipeline"
echo
echo "Syntax: FunOMIC -1 -2 -p -o -a -b -c -t"
echo "options:"
echo "1 path of paired end read1 (we recommend to apply kneaddata or other quality control first)"
echo "2 path of paired end read2 (we recommend to apply kneaddata or other quality control first)"
echo "p output prefix"
echo "o path of output directory"
echo "a path of UHGG bacterial database directory"
echo "b path of FunOMIC-T databases directory"
echo "c path of FunOMIC-P databases directory"
echo "t number of threads"
echo
}
######################################
########### SET VARIABLES ############
######################################
while getopts "1:2:p:o:a:b:c:t:h" flag
do
case "$flag" in
1) read1=$OPTARG;; # path of paired end read1 (we recommend to apply kneaddata or other quality control first)
2) read2=$OPTARG;; # path of paired end read2 (we recommend to apply kneaddata or other quality control first)
p) prefix=$OPTARG;; #output prefix
o) outdir=$OPTARG;; #path of output directory
a) bactdb=$OPTARG;; #path of UHGG bacterial database
b) taxadb=$OPTARG;; #path of FunOMIC-T databases
c) protdb=$OPTARG;; #path of FunOMIC-P databases
t) threads=$OPTARG;; #number of threads
h) Help
exit;; #display help
esac
done
# ##############################################
# ######### DECONTAMINATION WITH UHGG! #########
# ##############################################
printf "Start removing bacterial reads for ${prefix}\n"
if [[ ! -f $outdir ]]
then # When the pipeline has been executed properly and the script recognizes the output
mkdir $outdir
if [[ ! -f $outdir/tmp ]]
then mkdir $outdir/tmp
fi
if [[ ! -f $outdir/taxonomic_profiling ]]
then mkdir $outdir/taxonomic_profiling
fi
if [[ ! -f $outdir/functional_profiling ]]
then mkdir $outdir/functional_profiling
fi
fi &>/dev/null
bowtie2 -p $threads -x $bactdb/uhgg \
-1 $read1 \
-2 $read2 \
-S $outdir/tmp/${prefix}_Bact.sam 2> ${outdir}/bact_decontam.log
samtools view -b -f 4 $outdir/tmp/${prefix}_Bact.sam > $outdir/tmp/${prefix}_noBact.bam
samtools sort -n $outdir/tmp/${prefix}_noBact.bam -o $outdir/tmp/${prefix}_noBact_sorted.bam &>/dev/null
samtools fastq -@ 8 $outdir/tmp/${prefix}_noBact_sorted.bam \
-1 $outdir/${prefix}_noBact_1.fastq.gz \
-2 $outdir/${prefix}_noBact_2.fastq.gz \
-0 /dev/null -s /dev/null -n &>/dev/null
# ##############################################
# ########### TAXONOMICAL PROFILING! ###########
# ##############################################
printf "Start taxonomic annotation for ${prefix}\n"
bowtie2 -p $threads -x $taxadb/FunOMIC.T.v1 \
-1 $outdir/${prefix}_noBact_1.fastq.gz \
-2 $outdir/${prefix}_noBact_2.fastq.gz \
-S $outdir/tmp/${prefix}.sam 2> $outdir/taxonomic_profiling/log
# filter hits with q-score over 30 and coverage over 80
samtools view -Sq 30 $outdir/tmp/${prefix}.sam > $outdir/tmp/${prefix}.30.sam
samtools view -bSq 30 $outdir/tmp/${prefix}.sam > $outdir/tmp/${prefix}.30.bam
samtools sort -o $outdir/tmp/${prefix}.30.sorted.bam $outdir/tmp/${prefix}.30.bam &>/dev/null
coverageFilter.py -i $outdir/tmp/${prefix}.30.sam -o $outdir/tmp/${prefix}.30.c80.list &>/dev/null
samtools view $outdir/tmp/${prefix}.30.sorted.bam | grep -f $outdir/tmp/${prefix}.30.c80.list > $outdir/tmp/${prefix}.30.c80.sam
samtools view -bt $taxadb/FunOMIC.T.v1.fasta.fai -o $outdir/tmp/${prefix}.30.c80.bam $outdir/tmp/${prefix}.30.c80.sam &>/dev/null
samtools index $outdir/tmp/${prefix}.30.c80.bam &>/dev/null
samtools idxstats --threads $threads $outdir/tmp/${prefix}.30.c80.bam > $outdir/tmp/${prefix}.idxstats.txt
awk '$3 > 0' $outdir/tmp/${prefix}.idxstats.txt | sed 's/|/\t/g' | awk '{print $1,$4}' > $outdir/tmp/${prefix}"_assembly_hits.txt"
buglist.py -i $outdir/tmp/${prefix}"_assembly_hits.txt" -t $taxadb/taxonomy.txt -o $outdir"/taxonomic_profiling"/${prefix}"_buglist_stratified.txt" &>/dev/null
#### rm temporary files
rm -r $outdir/tmp/
##############################################
########### FUNCTIONAL PROFILING! ############
##############################################
printf "Starting functional annotation for ${prefix}\n"
cd ${outdir}/functional_profiling
mkdir $outdir/functional_profiling/tmp/
zcat $outdir/${prefix}_noBact* > $outdir/functional_profiling/joined.fastq
printf "merged clean reads stored in: $outdir/functional_profiling/joined.fastq \n"
diamond blastx -b5 -c1 --query-cover 95 \
--id 99 -d $protdb"/FunOMIC.P.v1.dmnd" \
-q $outdir"/functional_profiling"/joined.fastq --outfmt 6 \
-o $outdir/"functional_profiling"/${prefix}.func.out 2> $outdir/functional_profiling/log
# Filter alignment hits: coverage>80%; query length>100bp; sort by perc_identity, length, bit score, remove redundancy
awk '$4 >= 20' $outdir/"functional_profiling"/${prefix}.func.out | sort -k3,3nr -k4,4nr -k12,12nr | awk '!a[$1]++' > $outdir/"functional_profiling"/${prefix}.func.filtered.out
if [[ -f $outdir/"functional_profiling"/${prefix}.func.filtered.out ]]
then # When the pipeline has been executed properly and the script recognizes the output
printf "removing secondary files files for ${p}\n\n"
rm $outdir/functional_profiling/${prefix}.func.out
rm $outdir/functional_profiling/*.fastq
rm $outdir/functional_profiling/${prefix}.hist*
rm -r $outdir/functional_profiling/tmp
else echo "No output generated"
fi &>/dev/null
# functional_profiling.py -i $outdir/"functional_profiling"/${prefix}.func.filtered.out\
# -c $protdb/id_to_clean.txt -a $protdb/nju_id_protein_ec.txt \
# -o1 $outdir/"functional_profiling"/${prefix}"_protein.csv"\
# -o2 $outdir/"functional_profiling"/${prefix}"_EC.csv" &>/dev/null
get_pwy_abundance_v1.py -i $outdir/"functional_profiling"/${prefix}.func.filtered.out\
-s $protdb \
-o1 $outdir/"functional_profiling"/${prefix}"_pwy_abundance.csv"\
-o2 $outdir/"functional_profiling"/${prefix}"_pwyClass_abundance.csv"\
-o3 $outdir/"functional_profiling"/${prefix}"_pwyType_abundance.csv"\
-o4 $outdir/"functional_profiling"/${prefix}"_full_annotation.csv" #&>/dev/null