之前写过一个【基因id转换】的脚本,是针对fasta文件的。
一些软件对基因id的长度有一定限制,改完基因id跑完软件之后,发现要把结果文件中的id再改回来并不像改fasta文件那么简单。
下面举了一种结果文件的例子,其实还有很多种(如基因树文件中的id要改的情况...),关键就在于有一个原有id和改后id的对照表。
1_wrong_name.txt
、a_wrong_name.txt
:需要改基因id的文件
$ cat 1_wrong_name.txt
white gene1,gene2,gene3
black gene4,gene5
gray gene11,gene15
$ cat a_wrong_name.txt
red genea,geneb,genec
blue gened,genee
purple geneab,genedd
1_gene_id.txt
、a_gene_id.txt
:基因id对应表
$ cat 1_gene_id.txt
gene1 ID1
gene2 ID2
gene3 ID3
gene4 ID4
gene5 ID5
gene11 ID11
gene15 ID15
$ cat a_gene_id.txt
genea ida
geneb idb
genec idc
gened idd
genee ide
geneab idab
genedd iddd
for i in `echo 1 a`
do
j=`echo ${i}_wrong_name.txt`
awk '{print $2}' $j | sed 's/,/\n/g' > gene.txt
grep -wFf gene.txt ${i}_gene_id.txt | awk '{print $0,length($1)}' | sort -k 3nr | awk -va=$j 'NR==1{print "sed danyinhaos/"$1"/"$2"/gdanyinhao "a" | \\"}NR>1{print "sed -e danyinhaos/"$1"/"$2"/gdanyinhao | \\"}END{print "grep -v ^$ > changed_"a}' | sed "s/danyinhao/'/g" > change_gene.sh
bash change_gene.sh
rm gene.txt change_gene.sh
done
changed_1_wrong_name.txt
、changed_a_wrong_name.txt
:改完基因id之后的文件
$ cat changed_1_wrong_name.txt
white ID1,ID2,ID3
black ID4,ID5
gray ID11,ID15
$ cat changed_a_wrong_name.txt
red ida,idb,idc
blue idd,ide
purple idab,iddd
代码先提取结果文件基因id,再
grep
id对照表,确定要改的id,利用awk
生成一个含有sed
命令的脚本。
也可以不麻烦提取id,直接把所有基因id都放进脚本中,但几万个基因实在是...
网友评论