写在前面
##先mkdb好库
-rw-r--r-- 1 mkdb2020.sh
-rw-r--r-- 1 Paorum.fa
-rw-r--r-- 1 Paorum.fa.nhr
-rw-r--r-- 1 Paorum.fa.nin
-rw-r--r-- 1 Paorum.fa.nog
-rw-r--r-- 1 Paorum.fa.nsd
-rw-r--r-- 1 Paorum.fa.nsi
-rw-r--r-- 1 Paorum.fa.nsq
##再准备好脚本
-rw-r--r-- 1 work.sh
-rw-r--r-- 1 sp.lst
-rw-r--r-- 1 exonrate.sh
-rw-r--r-- 1 queries.fas
报错
print() on closed filehandle OUT at /zfssz3/NASCT_BACKUP/MS_PMO2017/lv/lv/software/ysoftware/bin/gene_predict/call_genewise.pl line 165, <IN> line 1.
print() on closed filehandle OUT at /zfssz3/NASCT_BACKUP/MS_PMO2017/lv/lv/software/ysoftware/bin/gene_predict/call_genewise.pl line 169, <IN> line 2.
- 报错与
sp.lst是否加.fa无关 - blast新旧版本不一致导致序列ID前加入了不可识别的
gb| |以致后续不能提取序列nuc，故blastfoamtdb和后面的都需保持一致版本
$ cat Aciculata.pos
L1998_T1/1_T_WB-D1 gb|AZMS01S00010088.1| + 70717 111986
L1998_T1/1_T_WB-D2 gb|AZMS01S00083152.1| + 54377 124604
L2916_T1/1_T_WB-D1 gb|AZMS01S00083152.1| + 54272 112779
- line 165不能print OUT 是因为
id这个变量没赋值成功导致的，进而推断是call_genewise.pl的输入文件有误
142 #################################################################################################
143 sub sequence_parser {
144 my $in_file = shift;
145 my $out_dir = shift;
146 my $type = shift;
147 open IN, $in_file;
148 open OUT, $in_file;
149 my $out_file;
150 while (<IN>) {
151 if (/^\s*>/) {
152 close OUT;
153 my $id;
154 if ($type eq "nuc") {
155 $id = (split /\s+/)[0];
156 $id =~ s/^>//;
157 $id = (split /##/, $id)[0];
158 } elsif ($type eq "pep") {
159 $id = (split /\s+/)[0];
160 $id =~ s/^>//;
161 }
162 $out_file = "$out_dir/$id.fa";
163 $out_file =~ s/\|/_/g;
164 open OUT, ">$out_file";
165 print OUT;
166 } else {
167 s/^\s+$//;
168 s/^\s+|\s+$//g;
169 print OUT uc($_), "\n" if $_;
170 }
171 }
172 close IN;
perl /zfssz3/NASCT_BACKUP/MS_PMO2017/lv/lv/software/ysoftware/bin/gene_predict/call_genewise.pl --pep $sp/$sp.pep --nuc $sp/$sp.nuc --list $sp/$sp.strandList --key $sp --out $sp/ --num $jobNum
cd ./$sp/result
- 检查各输入文件，发现各个ID上都带有/符号，推测是此原因，故而进行替换为_，后解决
L9283_T1/1_T_V-D4 -
L9283_T1/1_T_V-D5 -
Aciculata.strandList (END)
sed 's:/:xie:g' queries.fas > queries.fas3
sed 's:|:shu:g' queries.fas3 > queries.fas4
更新问题
最后预测得到的位置信息解释>DesHic_hcG_CAB89497.1-D6##HiC_scaffold_14##62993454##63035710(30257-12000)
Query: hcG_CAB89497.1-D6
Target: hcG_CAB89497.1-D6##HiC_scaffold_14##62993454##63035710 HiC_scaffold_14:62993454-63035710:- len=42257:[revcomp]
Model: protein2genome:local
Raw score: 2103
Query range: 2 -> 627
Target range: 30257 -> 12000
-
##62993454##63035710为原scaffold上的绝对碱基所在位置区间，是经过延伸截取出来做精细预测全长的
homologGenePredict.sh $pep.m8.solar.cor.idAdd.selected $pep $ref DesHic 12000 20 1
perl /zfssz3/NASCT_BACKUP/MS_PMO2017/bin/gene_predict/get_pos.pl $solar >$sp/$sp.pos
perl /zfssz3/NASCT_BACKUP/MS_PMO2017/bin/gene_predict/extract_sequence.pl --pos $sp/$sp.pos --fasta $genome --extent $extend >$sp/$sp.nuc
perl /zfssz3/NASCT_BACKUP/MS_PMO2017/bin/gene_predict/prepare_pep.pl $sp/$sp.pos $pep >$sp/$sp.pep
awk '{print $1 " "$3}' $sp/$sp.pos >$sp/$sp.strandList
perl /zfssz3/NASCT_BACKUP/MS_PMO2017//bin/gene_predict/call_genewise.pl --pep $sp/$sp.pep --nuc $sp/$sp.nuc --list $sp/$sp.strandList --key $sp --out $sp/ --num $jobNum
-
(30257-12000)为截取的基础上的确定位置，也就是##62993454在这个计算中是为第一位 -
请注意CDS DNA序列结果和核酸序列结果的不同
cds #预测得到去除了内含子的结果
final # 预测包含所有外显子和内含子信息
nuc # 预测前延伸截取的核酸目标区域，注意不是target的预测得出的核酸全长
pep # 预测前query的氨基酸序列
-
请注意
-是指反向互补序列编码，DNA编码原则-都是指反向互补编码，不只是反向！！！
Query range: 2 -> 627 Target range: 30257 -> 12000
less -i打开后是不区分大小写查找
CACGGTCATGGTCATCAT是GTGCCAGTACCAGTAGTA的反向互补序列
网友评论