josuebarrera
commented
9 months ago Looks very straightforward, thanks for adding this solution!
Closed Eulium closed 9 months ago
josuebarrera
commented
9 months ago Looks very straightforward, thanks for adding this solution!
Porpossale to fix bug in Erassignmet possible resulting in missing lineage information
BUG encountered:
When using NCBI's NR database with Diamond the output file can contain rows with multiple tax ids. These are separated by a ';' (see below, last column). In line 50 of the
Erassignmentscript the diamond output is piped into multiple awk and sed commands which extract the tax id column and use's its entries as key. In rows with multiple tax id's the ';' separated values will be used as a key in total, resulting in no matching rows from the phylostrata file. Any tax id that is only present in rows with multiple tax id's will not be used in later procedures of Erassignment and may cause wrong or incomplete results as phylostrat linage information matching tax id's found only in rows with ';' would not be present in thetmp_phylostratafile.Excerpt from a Diamond output file which will reproduce the bug:
Buggy code
Line 50
awk -v GENE="$GENE" '{ if ($1 == GENE) print $NF }' ${TMP_PATH}/tmp_${GENE}.bout | sort -u -T ${TMP_PATH} | sed '/^[[:space:]]*$/d' | awk -F "," 'FNR==NR{ a[$1]=$0;next } ($1 in a)' - ${NCBITAX} | grep "$OLDEST" > %{TMP_PATH}/tmp_phylostrataProposed change
Add
sed 's/;/\n/gcommand into the pipeline to split any ';' separated tax id's into multiple rows. This will result in them being split into single lines and properly used as single keys in the followingawkcommand Line 50awk -v GENE="$GENE" '{ if ($1 == GENE) print $NF }' ${TMP_PATH}/tmp_${GENE}.bout | sed 's/;/\n/g' | sort -u -T ${TMP_PATH} | sed '/^[[:space:]]*$/d' | awk -F "," 'FNR==NR{ a[$1]=$0;next } ($1 in a)' - ${NCBITAX} | grep "$OLDEST" > %{TMP_PATH}/tmp_phylostrata