需要解决的问题：

准确获取基因X的某个保守区域H的蛋白序列及核酸序列，查看该段序列蛋白层面及核酸层面的变化情况。

分析：

需要获取的基因X的某个保守区域序列H中有几个氨基酸非常保守，先进行序列比对，定位该段保守序列，提取该段序列，考虑到取到的序列中含有gap的情况，需要将获取的序列（去除“-”）重新在物种W的蛋白序列中定位再次提取，保证氨基酸数目的一致；定位得到的蛋白序列起始终止位置的坐标乘以3即得到核酸序列的坐标信息，从物种W的CDS序列文件中提取该保守区段H的核酸序列。

解决方案：

1、使用muscle将基因X的蛋白序列进行比对，并在Jalview中进行可视化

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

基因X序列比对

序列过长过短都会造成比对结果中产生gap，手动查看并删除造成很长gap的序列；
Ctrl+F定位该保守区段在比对文件中的位置，如上图确认PGRTDN出现的起始位置为123；
在Jalview中导出筛选之后的序列比对文件。
防止出错，去除比对文件的换行符

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比对文件中按保守区域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("请输入保守区域的起始位置：")
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
请输入保守区域的起始位置：123

输出文件为W_X_filter_aln_H.pep

3、从比对文件中提取的序列可能有gap，因此可能不足18个氨基酸，因此需要在原始蛋白序列文件中重新提取保证所有保守区域都是提取出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