使用 hmm文件搜索

hmmsearch --cut_tc --domtblout uboxhmmsearch.out U-box.hmm GCF_007655135.1_ASM765513v2_protein.faa

根据e值提取序列名称

grep -v "#" uboxhmmsearch.out | awk '$7 < 1e-20' | cut -f1 -d " " | sort -u > ubox_qua_id.txt

提取序列

seqtk subseq GCF_007655135.1_ASM765513v2_protein.faa ubox_qua_id.txt > ubox_qua_id.faa

序列比对后构建模型
hmmbuild NBS-ARC_qua.hmm NBS-ARC_qua.aln

使用构建的模型重新搜索，并根据e值提取序列名称

grep -v "#" ubox_secondsearch.out | awk '1}' | sort -u >ubox_second_id.txt

使用hmm模型搜索到序列
wc -l ../uboxid.txt
96 ../uboxid.txt
blast搜索到的序列

wc -l ../uboxblast/atuboxblastid.txt

1136 ../uboxblast/atuboxblastid.txt
合并以上id，获取序列

seqtk subseq ../GCF_007655135.1_ASM765513v2_protein.faa uboxsortid.txt > uboxsortid.faa

smart 检测序列结构域，获得含有ubox结构域的序列id。smart网站也可以检测pfam网站提供的结构域数据，smart 得到的结果有32条序列含有结构域。但是blast检测到的序列多达上千条，

E1 -3 分别为激活 藕联 连接酶 ，
将蛋白序列中 ubox domin like序列再pfam检测其结构域，有的序列并不含有ubox结构域信息
1.1 使用hmm文件搜索蛋白序列
中文文章不严谨，关于u-box的蛋白保守域pfam id 出错

# hmmsearch --cut_tc --domtblout uboxhmm.out U-box.hmm GCF_007655135.1_ASM765513v2_protein.faa 


# grep -v "#" uboxhmm.out | awk '{print $1}' >uboxhmmid.txt
搜索序列并提取序列id

1.2 blast 鉴定 序列

# cat 2kinduboxweb.txt  | sort | uniq -cu
1 AT2G40640
      1 AT3g47820
      1 AT5G05230
      1 AT5g57035

因为显示有2个ubox基因家族的网页，比较了2个网页的 locus信息，有四个locus 不是共有的。而且https://www.arabidopsis.org/browse/genefamily/plantubox.jsp 网页显示的基因家族只有63个
根据[拟南芥ubox基因家族](https://www.arabidopsis.org/browse/genefamily/plantubox.jsp) 提供的序列号
从UniProt[https://www.uniprot.org/](https://www.uniprot.org/)下载ubox的蛋白序列
对蛋白序列进行建库
# /root/app/blast+/ncbi-blast-2.10.0+/bin/makeblastdb -in GCF_007655135.1_ASM765513v2_protein.faa -dbtype prot -title GCF_007655135.1_ASM765513v2_protein -parse_seqids -logfile logfile.out

# /root/app/blast+/ncbi-blast-2.10.0+/bin/blastp  -query at63ubox_uniprot.faa -db GCF_007655135.1_ASM765513v2_protein_blast -out ubox.blast -outfmt 6 -evalue 1e-3

# wc -l ubox.blast
11041 ubox.blast

# cat ubox.blast | awk '{print $2}'|sort -u >blastp_ubox_id.txt
目测 hmm效果不大，直接使用blast就行，
base) [root@al hmmandblast]# wc -l ubox_all_id.txt 
1256 ubox_all_id.txt
(base) [root@al hmmandblast]# wc -l ../uboxhmm/uboxhmmid.txt
96 ../uboxhmm/uboxhmmid.txt
(base) [root@al hmmandblast]# wc -l ../uboxblast/blastp_ubox_id.txt
1160 ../uboxblast/blastp_ubox_id.txt
(base) [root@al hmmandblast]# cat ubox_all_id.txt |sort -u | wc -l
1160
(base) [root@al hmmandblast]# cat ubox_all_id.txt |sort -u >ubox_all_sort_id.txt 

# cat ubox_all_id.txt | sort | uniq -c |awk '$1==2'|wc -l
96  2中方式有96 个id 相同
# grep -v ">" ubox_all_sort_id.faa | awk -F "" '{print $1}'|sort -u
M
Q
S


>XP_031375845.1 第一个氨基酸 S
>XP_031390624.1 .......................Q
其他都是M，两条序列在smart no match

将结果为DOMAIN=Ubox的序列id提取为,id含有两个,这个基因有2个mrna variant，

XP_031385043.1
XP_031396537.1
XP_031393056.1
XP_031402432.1
XP_031401235.1
XP_031400178.1
XP_031404116.1
XP_031386047.1
XP_031395601.1
XP_031377406.1
XP_031394767.1
XP_031399446.1
XP_031391630.1
XP_031400912.1
XP_031398092.1XP_031398091.1
XP_031373647.1
XP_031405927.1
XP_031391114.1
XP_031397642.1
XP_031398133.1
XP_031379370.1
XP_031373608.1
XP_031402110.1
XP_031395361.1
XP_031389245.1
XP_031400164.1
XP_031403345.1
XP_031397910.1
XP_031391943.1XP_031391944.1
XP_031407253.1
XP_031389227.1

XP_031398092.1XP_031398091.1 的差异，转录为2种mrna 翻译出一种蛋白质

mrna variant 12

XP_031382446.1仅在pfam检测到u-box结构域，在smart检测到u-box结构域的序列都在pfam检测到结构域

明天用pfam检测所有序列结构域。慢，下载pfam数据库，数据库10个g下载失败，检测出的序列比文献报道 的少，再次用拟南芥和水稻的序列搜索，在http://rice.plantbiology.msu.edu/index.shtml 水稻数据库先搜索含有u box 基因的序列，这个网站太卡了，31条，locus含有染色体信息
使用拟南芥和水稻蛋白序列搜索，共搜索到1183条，但使用拟南芥的搜索到1160条

# cat at_os_ubox.blast |awk '{print $2}' |sort |uniq | wc -l
1183

单纯用拟南芥和 用（拟南芥与水稻）预测到相同的序列id

[root@al ubox]# cat ubox_all_id_from_at_and_atos.txt |sort|uniq -c |awk '$1==1'
      1 
      1 XP_031391943.1XP_031391944.1
      1 XP_031391944.1XP_031391943.1
      1 XP_031398091.1XP_031398092.1
      1 XP_031398092.1XP_031398091.1

pfam单独预测、拟南芥单独blast、拟南芥与pfam（重复id为pfam得到所有id）、拟南芥水稻加pfam 经smart可以检测到共同数量的目标结构域序列。
单独对hmmseach搜索到的序列提取低e值序列比对build重新搜索还没试。

# cat ubox_all_id_from_at_and_atos.txt |sort|uniq -c |awk '$1==3'|wc -l
33

根据低e序列重新使用hmmsearch搜索得到167条序列，但与拟南芥重合id 是140，有再次hmmsearch 得到27条序列没有通过blast得到，可能有因为这27条序列与拟南芥差别很大，第一次搜索的所有序列是否第二次全部搜索到

# grep -v "#" ubox_hmm_senond.out |awk '{print $1}' |sort -u|wc -l
167
# cat ubox_hmm_senond_sort.id /root/genefamily/ubox/uboxblast/at_ubox/blastp_ubox_id.txt |sort|uniq -c |awk '{print$1}'|sort|uniq -c
   1047 1
    140 2

# cat /root/genefamily/ubox/uboxhmm/uboxhmmid.txt|sort -u |wc -l
96
# cat /root/genefamily/ubox/uboxhmm/uboxhmmid.txt|wc -l
96
# cat ubox_hmm_senond_sort.id
/root/genefamily/ubox/uboxhmm/uboxhmmid.txt |sort|uniq -c |awk '{print$1}'|sort|uniq -c
     77 1
     93 2     所以第一次 搜索到的在第二次并没有全部搜索到。
可以用comm 在第二次搜索到的id 但是不在拟南芥blast到的id里，应该有27 条
查看一下这27条是否有结构域，，27条序列不含特定结构域不过包括DOMAIN=Pfam:zf-RING_UBOX

将33条序列提取出来，并重新在smart预测，花了将近一个下午，学会了将smart的结果提取出来,wuwuwu,加油，
这个是需要知道含有结构域的序列id，再次单独预测，然后提取，需要把除了id、domain的内容删除，不过直接针对smart的文本文件也可以，

import os 
os.getcwd()
'F:\\papers\\genefamily\\ubox\\at_os_hmm_id_sort'
os.chdir("F:/papers/genefamily/ubox/at_os_hmm_id_sort/")
fr1 = open ("at_os_hmm_ubox_id_sort_smart_id_faa_smart.txt","r")
fr2 = open("out.txt","w")
fr2.write("%s\t%s\t%s\t%s\t"%("accession_num","domain_type","start","end"))
fr2.write("\n")
1
smart = {}
for line in fr1:
    if line.strip().split("=")[0] == "USER_PROTEIN_ID":
        proids = line.strip().split("=")[1]
        smart[proids] = {}
    else:
        if line.strip().split("=")[0] == "DOMAIN":
            domainid = line.strip().split("=")[1]
            smart[proids][domainid] = {}
        else:
            k = line.strip().split("=")[0]
            v = line.strip().split("=")[1]
            smart[proids][domainid][k] = v
        
smart["XP_031396537.1"]["Ubox"]["EVALUE"]
'3.3627937278163e-10'
for a in smart:
    dt = "Ubox_smart"
    st =str(smart[a]["Ubox"]["START"])
    en = str(smart[a]["Ubox"]["END"])
    fr2.write("%s\t%s\t%s\t%s\t"%(a,dt,st,en))
    fr2.write("\n")

因为也可以使用很多物种的比对文件构建hmm,重新构建比对文件合并blast序列，搜索到相同数量的序列，
在SMART网站预测全蛋白组，但是只有四条序列还有ubox，共检测到10000多条序列。

20200320

比对参数

mege的比对文件导出后出错，应该再origin tree 里到处树文件。在bootstrap树导出出错，用iqtree构建可能比较准确

20200321
越看论文越感觉难，

从gff文件中得到 蛋白id与基因id的对应关系

# grep -v "#" GCF_007655135.1_ASM765513v2_genomic.gff |awk -F "[\t=;,:]" 'BEGIN{OFS="\t"}$3=="CDS"{print$15,$17}'|grep -f /root/genefamily/ubox/hmmandblast/at_os_hmm_ubox/at_os_hmm_ubox_id_sort_smart_id.txt |sort -u >proteinid_geneid.txt

同样提取id 的染色体及起始终止位置 正负链，及染色体长度等信息，将上述内容整合到一张表。

import os
os.getcwd()
'C:\\Users\\Acer'
os.chdir("C:/Users/Acer/Desktop/python/chormosome_location/")
proid_genid = open("proteinid_geneid.txt","r")
gene_info = open("33gene_info.txt","r")
chromo = open("8chromosome_info.txt","r")
out = open("out.txt","w")
proid_genid_dict = {}
gene_info_dict = {}
chromo_id_dict = {}
chromo_dict = {}
#整理蛋白id与基因id
for line in proid_genid:
    protein_id = line.strip().split("\t")[1]
    gene_id = line.strip().split("\t")[0]
    proid_genid_dict[protein_id] = gene_id
    
proid_genid_dict["XP_031373608.1"]
'116188407'
#基因的起始位置等信息信息
for line in gene_info: 
    gene_info_dict[line.strip().split("\t")[4]] = line.strip().split("\t")[0:4]
    chromo_id_dict[line.strip().split("\t")[4]] = line.strip().split("\t")[0]
    
gene_info_dict["116192860"]
['NC_045127.1', '52940497', '52944484', '-']
for line in chromo:
    chromo_dict[line.strip().split("\t")[0]] = line.strip().split("\t")[1:]
    
for a,b in  proid_genid_dict.items():
    out.write("%s\t%s\t"%(a,b))
    for gene_info_dict_value in gene_info_dict[b]:
        out.write("%s\t"%(gene_info_dict_value))
    chrom_id = chromo_id_dict[b]
    for chromo in chromo_dict[chrom_id]:
        out.write("%s\t"%(chromo))
    out.write("\n")

    
out.close()

根据基因的起始位置作图需要考虑正负链么，基因在负链上时，gff文件里的起始位置就时基因在负链上基因的结束位置

基因在负链上

不同染色体上的数量

# grep -v "#" gene_info.txt |awk '{print $9}'|sort|uniq -c
      1 chromosome=1
      1 chromosome=2
      4 chromosome=3
      6 chromosome=4
      8 chromosome=5
      7 chromosome=6
      4 chromosome=7
      2 chromosome=8

gff文件中，基因id 转录本id 蛋白 id 不是加一
号，
从gtf文件根据蛋白id提取前8列内容及把 transcript_id 更换为基因的简写
突然想到要是直接把某行的内容分割出来，只要考虑字符串两端的符号就行了，啊啊啊，不需要把每一行的内容尽可能地分割出来

一个下午，勉强可以从gtf文件中获取需要的序列，并更改转录本名称为自定义

import os
import re
os.chdir("C:/Users/Acer/Desktop/python/gene_struc/")
gtf = open("GCF_007655135.1_ASM765513v2_genomic.gtf","r")
gene_info = opimport osen("gene_info.txt","r")  #蛋白与mrna，基因id对应文件
out = open("out.txt","w")
gene_info_dict = {}
for line in gene_info:
    gene_info_dict[line.strip().split("\t")[2]] = line.strip().split("\t")[10]   #以mrna id为键，自定义名称为值
for line in gtf:
    if line[0] == "#" :
        continue
    elif line.strip().split("\t")[2] == "gene":
        continue
    transcript_id  = line.strip().split("\"")[3]
    for i in gene_info_dict.keys():
        if(i == transcript_id):
            for aa in line.strip().split("\"")[0:3]: #transcript_id 以"为分割
                print(aa)
                out.write("%s\""%(aa))
            out.write("%s\""%(gene_info_dict[transcript_id]))
            for bb  in line.strip().split("\"")[4:]:
                out.write("%s\""%(bb))
            out.write("\n")
print("good job")
out.close()

但是得到的基因结构图没有utr的内容，推测可能时提取的时候没有提取第三列为 gene 的那一行。

gtf含有第三列gene的行时

但不含有基因那一行，也有utr

不含有基因那一行，也有utr

可能是因为基因与外显子起始位置相同，而cds与基因外显子不同，才有utr，基因优势基因做标尺，没有基因是外显子做标尺。

测试数据 单独画exon gene 与exon 基因与cds