Commit ebe03bf2 authored by peguerin's avatar peguerin
Browse files

very efficient script to query UNIPROT id from DNA sequence

parent ad0813a9
colonne_rank_query () {
uniprot_blastx_out=$1
j="rien"
count=0
for i in `grep -ve "^#" ${uniprot_blastx_out} | cut -d $'\t' -f 1`;
do if [ "$i" == "$j" ];
then count=$((count+1));
else count=1 ;
fi ;
echo $count ;
j="$i";
done
}
for orgn in `echo diplodus mullus`;
do for typeanot in `echo CDS` ;
do blastx -db swissprot_nucl -query ${orgn}_${typeanot}_flanking.fasta -task blastx -outfmt "7 delim=, qacc sacc evalue bitscore qcovus pident" -max_target_seqs 5 > ${orgn}_${typeanot}_uniprot_blastx.out;
grep -ve "^#" ${orgn}_${typeanot}_uniprot_blastx.out | tr ":" "\t" | cut -d $'\t' -f 1,2,3,4,5,7 > ${orgn}_${typeanot}_uniprot_blastx.out.tsv;
colonne_rank_query ${orgn}_${typeanot}_uniprot_blastx.out > ${orgn}_${typeanot}_uniprot_blastx.out.colon_id;
paste ${orgn}_${typeanot}_uniprot_blastx.out.tsv ${orgn}_${typeanot}_uniprot_blastx.out.colon_id > ${orgn}_${typeanot}_uniprot_blastx_no_header.tsv;
cat <(echo -e "chrom\tpos\tuniprot\tEvalue\tbitscore\tidentity\trank") ${orgn}_${typeanot}_uniprot_blastx_no_header.tsv > ${orgn}_${typeanot}_uniprot_blastx.tsv;
rm ${orgn}_${typeanot}_uniprot_blastx.out.tsv ${orgn}_${typeanot}_uniprot_blastx.out.colon_id ${orgn}_${typeanot}_uniprot_blastx_no_header.tsv;
done
done
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment