需要解決的問題：

準確獲取基因X的某個保守區(qū)域H的蛋白序列及核酸序列雾叭，查看該段序列蛋白層面及核酸層面的變化情況落蝙。

分析：

需要獲取的基因X的某個保守區(qū)域序列H中有幾個氨基酸非常保守暂幼，先進行序列比對，定位該段保守序列旺嬉，提取該段序列，考慮到取到的序列中含有g(shù)ap的情況捐顷，需要將獲取的序列（去除“-”）重新在物種W的蛋白序列中定位再次提取雨效，保證氨基酸數(shù)目的一致；定位得到的蛋白序列起始終止位置的坐標(biāo)乘以3即得到核酸序列的坐標(biāo)信息徽龟，從物種W的CDS序列文件中提取該保守區(qū)段H的核酸序列。

解決方案：

1传透、使用muscle將基因X的蛋白序列進行比對极颓，并在Jalview中進行可視化

muscle -in W-X.fas -out W-X.fas.aln -maxiters 1000

基因X序列比對

序列過長過短都會造成比對結(jié)果中產(chǎn)生gap，手動查看并刪除造成很長gap的序列兵琳；
Ctrl+F定位該保守區(qū)段在比對文件中的位置，如上圖確認PGRTDN出現(xiàn)的起始位置為123闰围；
在Jalview中導(dǎo)出篩選之后的序列比對文件。
防止出錯碧查，去除比對文件的換行符

awk '!/^>/ { printf "%s", $0; n = "\n" } /^>/ { print n $0; n = "" } END { printf "%s", n } ' W_X_filter_aln.fas.aln > W_X_filter_aln.fas.aln.nochangeline

2校仑、寫腳本從aln比對文件中按保守區(qū)域H的比對起始終止位置提取序列并去除gap（18個氨基酸）

#!/usr/bin/env python
# -*- coding:utf-8 -*-
import sys
import re
File = sys.argv[1]
input_file = open(File,'r')
start = input("請輸入保守區(qū)域的起始位置：")
end = int(start) + 18
prefix = sys.argv[1].split(".")[0]
output_file = open ('%s_H.pep'%(prefix), 'w')
for line in input_file:
    line = line.rstrip()
    if not line: break
    if '>' in line:
        print(line,file = output_file)
    else:
        line = line[int(start):int(end)]
        if "-" in line:
            line = line.replace("-", "")
            print(line,file = output_file)
        else:
            print(line,file = output_file)
input_file.close()
output_file.close()

#運行命令
python extract_H_pepseq.py W_X_filter_aln.fas.aln.nochangeline
請輸入保守區(qū)域的起始位置：123

輸出文件為W_X_filter_aln_H.pep

3、從比對文件中提取的序列可能有g(shù)ap迄沫，因此可能不足18個氨基酸，因此需要在原始蛋白序列文件中重新提取保證所有保守區(qū)域都是提取出18個氨基酸泰佳；
用腳本確定這些保守序列在原始蛋白序列中的起始終止位置

#!/usr/bin/env python
# -*- coding:utf-8 -*-
import sys
import re
File1 = sys.argv[1]
File2 = sys.argv[2]
input_file1 = open(File1,'r')
input_file2 = open(File2,'r')
prefix = sys.argv[2].split(".")[0]
output_file1 = open ('%s.peppos'%(prefix), 'w')
output_file2 = open ('%s.CDSpos'%(prefix), 'w')
Dict1 = dict()
Dict2 = dict()
for line1 in input_file1:
    line1 = line1.rstrip()
    if not line1: break
    if '>' in line1:
        idline1 = line1
    else:
        Dict1[idline1] = line1
for line2 in input_file2:
    line2 = line2.rstrip()
    if not line2:break  
    if '>' in line2:
        idline2 = line2
    else:
        Dict2[idline2] = line2
for key2,value2 in Dict2.items():
    for key1,value1 in Dict1.items():
        if key2 == key1:
            start = int(value1.find(value2)) + 1
            end = int(start) + 17
            print(key2.replace(">",""),start,end,file = output_file1)
            print(key2.replace(">",""),int(start)*3,int(end)*3+3,file = output_file2)
input_file1.close()
input_file2.close()
output_file1.close()
output_file2.close()

#運行命令
python find_Hpep_pos.py ../W.pep W_X_filter_aln_H.pep

輸出文件為W_X_filter_aln_H.peppos和W_X_filter_aln_H.CDSpos
輸出文件的格式如下：
<ID號> <起始位置> <終止位置>

** 4逝她、以上準備工作做完后用TBtools提取序列并出seqlogo圖**
用TBtools提取序列

java -cp TBtools_JRE1.6.jar biocjava.bioIO.FastX.FastaIndex.SeqExtracterAccordingFAindex --inFA ../W.pep --regionList W_X_filter_aln_H.peppos --outFA W_X_H.pep
java -cp TBtools_JRE1.6.jar biocjava.bioIO.FastX.FastaIndex.SeqExtracterAccordingFAindex --inFA ../W.CDS --regionList W_X_filter_aln_H.CDSpos --outFA W_X_H.CDS

用TBtools出seqlogo圖

java -cp TBtools_JRE1.6.jar biocjava.bioDoer.seqLogo.makeSeqLogo --inFile W_X_H.pep --borderSize 5 --borderColor 0,0,0 --scaleIC true --OutGraph W_X_H.pep.png
java -cp TBtools_JRE1.6.jar biocjava.bioDoer.seqLogo.makeSeqLogo --inFile W_X_H.CDS --borderSize 5 --borderColor 0,0,0 --scaleIC true --OutGraph W_X_H.CDS.png

完成２谴贰！擒贸！

H蛋白序列

H核酸序列

Usage

java -cp TBtools_JRE1.6.jar biocjava.bioDoer.seqLogo.makeSeqLogo --inFile W_X_H.pep --borderSize 5 --borderColor 0,0,0 --OutGraph W_X_H.pep1.png

--scaleIC false

選擇 --onlyBorder的變化如下：

java -cp TBtools_JRE1.6.jar biocjava.bioDoer.seqLogo.makeSeqLogo --inFile W_X_H.pep --borderSize 5 --borderColor 0,0,0 --scaleIC true --onlyBorder true --OutGraph W_X_H.pep2.png

--onlyBorder true

改變x軸和y軸的間距

java -cp TBtools_JRE1.6.jar biocjava.bioDoer.seqLogo.makeSeqLogo --inFile W_X_H.pep --borderSize 5 --borderColor 0,0,0 --scaleIC true --xInterval 5 --yInterval 5 --OutGraph W_X_H.pep3.png

間距1--5