BLAST比对得到目标序列

BLAST比对得到目标序列

这里我们不用服务器,同样我们利用的是Git软件,同样可以进行blast比对

此处以拟南芥为例:

在拟南芥基因组文件夹中右键bash here

第一步:建库


smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ ls
Arabidopsis_thaliana.TAIR10.42.gff3  Arabidopsis_thaliana.TAIR10.cds.all.fa  Arabidopsis_thaliana.TAIR10.dna.toplevel.fa  Arabidopsis_thaliana.TAIR10.pep.all.fa  拟南芥基因组/

smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ mv Arabidopsis_thaliana.TAIR10.pep.all.fa all.pep.fa

smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ ls
all.pep.fa  Arabidopsis_thaliana.TAIR10.42.gff3  Arabidopsis_thaliana.TAIR10.cds.all.fa  Arabidopsis_thaliana.TAIR10.dna.toplevel.fa  拟南芥基因组/

smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ #建库,此处是以蛋白建库,并以蛋白进行比对

smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ makeblastdb -in all.pep.fa -dbtype prot -title all.pep.fa


Building a new DB, current time: 04/10/2019 10:45:21
New DB name:   E:\genome_wide\?????\all.pep.fa
New DB title:  all.pep.fa
Sequence type: Protein
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 48321 sequences in 4.14913 seconds.

第二步:比对

此处是在Uniprot(https://www.uniprot.org/)找的一条wrky序列。





$ ls##建库后就会出现*.phr;*.pin;*.psq这三个文件
all.pep.fa      all.pep.fa.pin  Arabidopsis_thaliana.TAIR10.42.gff3     Arabidopsis_thaliana.TAIR10.dna.toplevel.fa  拟南芥基因组/
all.pep.fa.phr  all.pep.fa.psq  Arabidopsis_thaliana.TAIR10.cds.all.fa  wrky.fa
#设置E值为1e-10
smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ blastall -i wrky.fa -d all.pep.fa -p blastp -e 1e-10 -m 8 -o nbs_out_blast.txt
smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ ls
all.pep.fa      all.pep.fa.pin  Arabidopsis_thaliana.TAIR10.42.gff3     Arabidopsis_thaliana.TAIR10.dna.toplevel.fa  wrky.fa
all.pep.fa.phr  all.pep.fa.psq  Arabidopsis_thaliana.TAIR10.cds.all.fa  nbs_out_blast.txt                            拟南芥基因组/

smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ less nbs_out_blast.txt
sp|Q6QHD1|WRK71_ORYSJ   AT1G80840.1     36.54   353     166     8       1       347     1       301     7.2e-046        181.8
sp|Q6QHD1|WRK71_ORYSJ   AT4G31800.1     34.80   342     177     8       11      346     9       310     2.4e-041        166.8
sp|Q6QHD1|WRK71_ORYSJ   AT2G25000.1     36.96   303     124     6       44      346     36      271     1.1e-038        157.9
sp|Q6QHD1|WRK71_ORYSJ   AT2G25000.3     36.88   301     123     6       46      346     17      250     3.3e-038        156.4
sp|Q6QHD1|WRK71_ORYSJ   AT2G25000.4     37.07   294     118     6       53      346     2       228     7.2e-038        155.2
sp|Q6QHD1|WRK71_ORYSJ   AT2G25000.2     37.07   294     118     6       53      346     2       228     7.2e-038        155.2
sp|Q6QHD1|WRK71_ORYSJ   AT4G31800.2     48.88   178     64      3       173     346     91      245     1.2e-037        154.5
sp|Q6QHD1|WRK71_ORYSJ   AT4G31800.3     34.23   222     127     5       11      230     9       213     4.1e-025        112.8
sp|Q6QHD1|WRK71_ORYSJ   AT4G01720.1     32.73   220     145     2       45      264     94      310     1.6e-024        110.9
sp|Q6QHD1|WRK71_ORYSJ   AT4G04450.1     61.97   71      27      0       191     261     290     360     6.6e-023        105.5
sp|Q6QHD1|WRK71_ORYSJ   AT4G04450.2     61.97   71      27      0       191     261     196     266     6.6e-023        105.5
sp|Q6QHD1|WRK71_ORYSJ   AT1G62300.1     28.89   225     148     2       49      261     156     380     1.1e-022        104.8
sp|Q6QHD1|WRK71_ORYSJ   AT4G22070.4     65.63   64      22      0       191     254     295     358     3.3e-022        103.2
sp|Q6QHD1|WRK71_ORYSJ   AT4G22070.3     65.63   64      22      0       191     254     295     358     3.3e-022        103.2
sp|Q6QHD1|WRK71_ORYSJ   AT4G22070.2     65.63   64      22      0       191     254     295     358     3.3e-022        103.2
sp|Q6QHD1|WRK71_ORYSJ   AT4G22070.1     65.63   64      22      0       191     254     295     358     3.3e-022        103.2
sp|Q6QHD1|WRK71_ORYSJ   AT1G18860.1     55.00   80      33      1       192     271     238     314     2.1e-021        100.5
sp|Q6QHD1|WRK71_ORYSJ   AT1G68150.1     59.21   76      28      1       192     267     234     306     2.8e-021        100.1
sp|Q6QHD1|WRK71_ORYSJ   AT5G15130.2     64.52   62      22      0       192     253     184     245     3.6e-021        99.75
sp|Q6QHD1|WRK71_ORYSJ   AT5G15130.1     64.52   62      22      0       192     253     226     287     3.6e-021        99.75
sp|Q6QHD1|WRK71_ORYSJ   AT2G04880.2     53.85   78      33      2       177     254     269     343     2.0e-019        93.97
sp|Q6QHD1|WRK71_ORYSJ   AT2G04880.2     47.62   63      31      2       191     253     109     169     7.5e-011        65.47
sp|Q6QHD1|WRK71_ORYSJ   AT2G04880.1     53.85   78      33      2       177     254     293     367     2.0e-019        93.97
sp|Q6QHD1|WRK71_ORYSJ   AT2G04880.1     47.62   63      31      2       191     253     109     169     7.5e-011        65.47
sp|Q6QHD1|WRK71_ORYSJ   AT2G23320.1     41.12   107     60      1       158     261     202     308     2.4e-017        87.04
sp|Q6QHD1|WRK71_ORYSJ   AT5G56270.1     48.19   83      41      2       191     272     485     566     2.4e-017        87.04
sp|Q6QHD1|WRK71_ORYSJ   AT5G56270.1     51.43   70      32      2       193     262     273     340     7.3e-014        75.48
sp|Q6QHD1|WRK71_ORYSJ   AT2G37260.1     53.33   75      33      2       194     268     166     238     4.1e-017        86.27
sp|Q6QHD1|WRK71_ORYSJ   AT2G37260.1     39.18   97      54      3       157     250     311     405     1.4e-012        71.25
sp|Q6QHD1|WRK71_ORYSJ   AT2G37260.2     53.33   75      33      2       194     268     84      156     4.1e-017        86.27
sp|Q6QHD1|WRK71_ORYSJ   AT2G37260.2     39.18   97      54      3       157     250     229     323     1.4e-012        71.25
sp|Q6QHD1|WRK71_ORYSJ   AT2G37260.3     53.33   75      33      2       194     268     86      158     4.1e-017        86.27
sp|Q6QHD1|WRK71_ORYSJ   AT2G37260.3     39.18   97      54      3       157     250     231     325     1.4e-012        71.25
sp|Q6QHD1|WRK71_ORYSJ   AT5G41570.1     47.19   89      44      2       162     250     69      154     4.1e-017        86.27
sp|Q6QHD1|WRK71_ORYSJ   AT1G69810.2     46.88   96      48      2       170     264     124     217     7.0e-017        85.50
sp|Q6QHD1|WRK71_ORYSJ   AT1G69810.1     46.88   96      48      2       170     264     182     275     7.0e-017        85.50
sp|Q6QHD1|WRK71_ORYSJ   AT4G18170.1     54.84   62      27      1       192     253     171     231     2.7e-016        83.57
sp|Q6QHD1|WRK71_ORYSJ   AT2G38470.1     43.43   99      47      4       155     253     332     421     3.5e-016        83.19
sp|Q6QHD1|WRK71_ORYSJ   AT2G38470.1     48.61   72      35      2       193     264     184     253     9.5e-014        75.10
sp|Q6QHD1|WRK71_ORYSJ   AT1G13960.1     52.11   71      33      1       191     261     407     476     3.5e-016        83.19
sp|Q6QHD1|WRK71_ORYSJ   AT1G13960.1     31.82   132     85      3       185     316     221     347     9.5e-014        75.10
sp|Q6QHD1|WRK71_ORYSJ   AT1G13960.2     52.11   71      33      1       191     261     254     323     3.5e-016        83.19
sp|Q6QHD1|WRK71_ORYSJ   AT1G13960.2     31.82   132     85      3       185     316     68      194     9.5e-014        75.10
sp|Q6QHD1|WRK71_ORYSJ   AT1G55600.1     50.63   79      34      2       194     272     308     381     3.5e-016        83.19
sp|Q6QHD1|WRK71_ORYSJ   AT3G01080.2     46.25   80      42      1       191     270     342     420     6.0e-016        82.42
sp|Q6QHD1|WRK71_ORYSJ   AT3G01080.2     46.67   75      30      3       194     268     206     270     3.4e-011        66.63
sp|Q6QHD1|WRK71_ORYSJ   AT3G01080.1     46.25   80      42      1       191     270     304     382     6.0e-016        82.42
sp|Q6QHD1|WRK71_ORYSJ   AT3G01080.1     46.67   75      30      3       194     268     168     232     3.4e-011        66.63
#输出ID
smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ awk '{print$2}' nbs_out_blast.txt > id.txt

有一个桌面版本,这个也是桌面版本但是不是hmmer(参考:如何利用HMMER鉴定基因家族成员),而是blast,通过建库得到相似的序列。如果你想用hmmsearch也是可以的哦:





序列比对?

muscle官网(https://www.drive5.com/)

详情可参考:小工具使用<今天工具很散很乱>

同样可以,仅仅需要添加到环境既可:


smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ hmmsearch -h
# hmmsearch :: search profile(s) against a sequence database
# HMMER 3.0 (March 2010); http://hmmer.org/
# Copyright (C) 2010 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: hmmsearch [options] <query hmmfile> <target seqfile>

where basic options are:
  -h : show brief help on version and usage

options directing output:
  -o <f>          : direct output to file <f>, not stdout
  -A <f>          : save multiple alignment of all hits to file <s>
  --tblout <f>    : save parseable table of per-sequence hits to file <s>
  --domtblout <f> : save parseable table of per-domain hits to file <s>
  --acc           : prefer accessions over names in output
  --noali         : don't output alignments, so output is smaller
  --notextw       : unlimit ASCII text output line width
  --textw <n>     : set max width of ASCII text output lines  [120]  (n>=120)

options controlling reporting thresholds:
  -E <x>     : report sequences <= this E-value threshold in output  [10.0]  (x>0)
  -T <x>     : report sequences >= this score threshold in output
  --domE <x> : report domains <= this E-value threshold in output  [10.0]  (x>0)
  --domT <x> : report domains >= this score cutoff in output

options controlling inclusion (significance) thresholds:
  --incE <x>    : consider sequences <= this E-value threshold as significant
  --incT <x>    : consider sequences >= this score threshold as significant
  --incdomE <x> : consider domains <= this E-value threshold as significant
  --incdomT <x> : consider domains >= this score threshold as significant

options controlling model-specific thresholding:
  --cut_ga : use profile's GA gathering cutoffs to set all thresholding
  --cut_nc : use profile's NC noise cutoffs to set all thresholding
  --cut_tc : use profile's TC trusted cutoffs to set all thresholding

options controlling acceleration heuristics:
  --max    : Turn all heuristic filters off (less speed, more power)
  --F1 <x> : Stage 1 (MSV) threshold: promote hits w/ P <= F1  [0.02]
  --F2 <x> : Stage 2 (Vit) threshold: promote hits w/ P <= F2  [1e-3]
  --F3 <x> : Stage 3 (Fwd) threshold: promote hits w/ P <= F3  [1e-5]
  --nobias : turn off composition bias filter

other expert options:
  --nonull2     : turn off biased composition score corrections
  -Z <x>        : set # of comparisons done, for E-value calculation
  --domZ <x>    : set # of significant seqs, for domain E-value calculation
  --seed <n>    : set RNG seed to <n> (if 0: one-time arbitrary seed)  [42]
  --tformat <s> : assert target <seqfile> is in format <s>>: no autodetection
  --cpu <n>     : number of parallel CPU workers to use for multithreads
#利用模型搜索目标序列
smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$hmmsearch --domtblout outfile.txt --cut_tc *.hmm all.pep.fa





smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ muscle -h
Invalid command line option "h"

MUSCLE v3.8.31 by Robert C. Edgar

http://www.drive5.com/muscle
This software is donated to the public domain.
Please cite: Edgar, R.C. Nucleic Acids Res 32(5), 1792-97.


Basic usage

    muscle -in <inputfile> -out <outputfile>

Common options (for a complete list please see the User Guide):

    -in <inputfile>    Input file in FASTA format (default stdin)
    -out <outputfile>  Output alignment in FASTA format (default stdout)
    -diags             Find diagonals (faster for similar sequences)
    -maxiters <n>      Maximum number of iterations (integer, default 16)
    -maxhours <h>      Maximum time to iterate in hours (default no limit)
    -html              Write output in HTML format (default FASTA)
    -msf               Write output in GCG MSF format (default FASTA)
    -clw               Write output in CLUSTALW format (default FASTA)
    -clwstrict         As -clw, with 'CLUSTAL W (1.81)' header
    -log[a] <logfile>  Log to file (append if -loga, overwrite if -log)
    -quiet             Do not write progress messages to stderr
    -version           Display version information and exit

Without refinement (very fast, avg accuracy similar to T-Coffee): -maxiters 2
Fastest possible (amino acids): -maxiters 1 -diags -sv -distance1 kbit20_3
Fastest possible (nucleotides): -maxiters 1 -diags
smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$  muscle -in test.fa -out test.out.fa

MUSCLE v3.8.31 by Robert C. Edgar

http://www.drive5.com/muscle
This software is donated to the public domain.
Please cite: Edgar, R.C. Nucleic Acids Res 32(5), 1792-97.

test 6 seqs, max length 924, avg  length 318
00:00:00      4 MB(0%)  Iter   1  100.00%  K-mer dist pass 1
00:00:00      4 MB(0%)  Iter   1  100.00%  K-mer dist pass 2
00:00:00      7 MB(0%)  Iter   1  100.00%  Align node
00:00:00      7 MB(0%)  Iter   1  100.00%  Root alignment
00:00:00      7 MB(0%)  Iter   2  100.00%  Root alignment
00:00:00      8 MB(0%)  Iter   3  100.00%  Refine biparts
00:00:00      8 MB(0%)  Iter   4  100.00%  Refine biparts
00:00:00      8 MB(0%)  Iter   5  100.00%  Refine biparts
00:00:00      8 MB(0%)  Iter   5  100.00%  Refine biparts
00:00:00      8 MB(0%)  Iter   6  100.00%  Refine biparts

smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ ls
all.pep.fa      all.pep.fa.pin  Arabidopsis_thaliana.TAIR10.42.gff3     Arabidopsis_thaliana.TAIR10.dna.toplevel.fa  nbs_out_blast.txt  test.out.fa  拟南芥基因组/
all.pep.fa.phr  all.pep.fa.psq  Arabidopsis_thaliana.TAIR10.cds.all.fa  id.txt                                       test.fa            wrky.fa

smyang2018@smyang2018-PC MINGW64 /e/genome_wide/拟南芥
$ less test.out.fa

相关文章:

如何利用HMMER鉴定基因家族成员

顺式作用元件可视化

iTOL绘制进化树

MCScanX-transposed简介、安装及使用

染色体定位

基因结构可视化之—我来也!

收集 | 序列提取工具

延伸学习:

收集 | 序列提取工具

ascp咋啦?如何从NCBI上下载sra数据?

服务器 | 环境变量配置

服务器上安装最新版R和Rstudio(2019/3/20)

你在用xshell,putty登陆?推荐一个小工具(Git)登陆

怎样使用浏览器进行搜索

小工具使用<今天工具很散很乱>

服务器 | 又又买了一个服务器

如若侵犯您版权或者您有疑惑,请及时联系我,谢谢!

smyang2018

跳至工具栏