Hi all,
I'm working with snakemake to make automated workflow on my analysis, and I'm stuck at BLAST stage..
What I want to do is BLASTp against pre-made database --> Extract blasted queries having above 70 % identity & header section --> combine them as a single output.
There is no problem on "rule RAWBLAST_Hs" stage, and I find out that it works perfectly with other samples if there are any queries above 70% identity. But it makes an error if there is no queries with above 70% identity (which means it makes a 0 sized file) and try to use it at the rule FLT_BLAST_OUTCOME_Hs. "rule RAWBLAST_Hs" also works well with this sample. I run manually on bash, but it doesn't make any error... I guess "awk" command makes an error somehow with 0 sized temp_file_Hs and try to use it, but I have no idea what's wrong with it and how to fix it.
Below is my snakemake command:
rule RAWBLAST_Hs:
input:
config["QUERY_PATH"]+"{filename}.faa"
output:
config["OUTPUT_PATH"]+"BLAST_outcome/{filename}_Hs.RAW"
threads:
30
resources:
nodes=1,
ppn=30,
mem_gb=20,
walltime=86400
shell:
"""
module load tools perl/5.30.2 ncbi-blast/2.12.0+
blastp {config[Hs_BLASTParams]} \
-db {config[HsapiensDataBaseProt]} \
-query {input} \
-out {output} \
-outfmt 7 \
-num_threads {threads}
"""
rule FLT_BLAST_OUTCOME_Hs:
input:
rules.RAWBLAST_Hs.output
output:
config["OUTPUT_PATH"]+"BLAST_outcome/{filename}_Hs.FILTERED"
params:
temp_file_Hs="temp_file_Hs.txt"
shell:
"""
awk '$3>70' {input} > {params.temp_file_Hs}
grep "Fields:" {input} | uniq > header.txt
if [ -s {params.temp_file_Hs} ]
then
cat header.txt {params.temp_file_Hs} > {output}
else
cat header.txt > {output}
fi
"""
Hi, thank you for your comment. I set the BLASTp outfmt as 7, so even if there is no matched queries, there are always headers (which means, there is always return value from grep command). But thank you again for your opinion again, cuz I did not know grep returns non-zero exit if there is no matches :)
Ok, see edited answer...
Okay, so I follow your edited answer:
awk '$3>70' Sample_1 | grep -v '#' > test.txt
test.txt is the file which exists but is empty, and if I do as your edited answer with one change:
if [ ! -s test.txt]; then echo $?; fi
I got "0", not "1"