通过简单数据熟悉Linux下生物信息学各种操作1

原地址
几点说明
1.非简单翻译，所有代码均可运行，为了辅助理解，基本每步代码都有结果，需要比较的进行了整合
2.原文中的软件都下载最新版本
3.原文中有少量代码是错误的，这里进行了修正
4.对于需要的一些知识背景，在这里进行了注释或链接到他人博客

一共4部分
通过简单数据熟悉Linux下生物信息学各种操作1
通过简单数据熟悉Linux下生物信息学各种操作2
通过简单数据熟悉Linux下生物信息学各种操作3
通过简单数据熟悉Linux下生物信息学各种操作4

1下载酵母基因组gff格式文件

wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.gff.gz

2安装sratoolkit

curl https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.6-1/sratoolkit.2.9.6-1-mac64.tar.gz
tar xzvf sratoolkit.2.9.6-1-mac64.tar.gz
echo export PATH=$PATH:~/src/sratoolkit.2.3.5-2-mac64/bin >> ~/.profile
source ~/.profile

3下载数据

prefetch SRR1553610
find ~/ncbi

自动下载到家目录下的ncbi文件夹

/Users/ucco/ncbi
/Users/ucco/ncbi/public
/Users/ucco/ncbi/public/sra
/Users/ucco/ncbi/public/sra/SRR1553610.sra

sra格式转变为fq格式

fastq-dump --split-files SRR1553610.sra

wc -l *.fastq
  879348 SRR1553610_1.fastq
  879348 SRR1553610_2.fastq
 1758696 total

一次下载多个文件

$ echo SRR1553608 > sra.ids
$ echo SRR1553605 >> sra.ids
$ prefetch --option-file sra.ids

其他几种下载方式，看

从ncbi下载sra数据的几种种方式

4 通过EDirect获取序列

4.1根据locus获取序列

efetch -db nucleotide -id KM233090 -format fasta > KM233090.fa

4.2 根据accession number获取序列

efetch -db nucleotide -id 667853062 -format fasta > 667853062.fa

ucco:yeast ucco$ cat KM233090.fa |wc
     270     277   19169
ucco:yeast ucco$ cat 667853062.fa |wc
     270     277   19169

4.3同时获取多个序列

efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta > multi.fa

看下header files

$ cat multi.fa |grep ">"
>KM233090.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3816, complete genome
>KM233066.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3769.2, complete genome
>KM233113.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3856.1, complete genome

4.4只获取其中一部分

efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta -seq_start 1 -seq_stop 10 > multi.fa

$ cat multi.fa 
>KM233090.1:1-10 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3816, complete genome
AAATTGTTAC
>KM233066.1:1-10 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3769.2, complete genome
GAATAACTAT
>KM233113.1:1-10 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3856.1, complete genome
CGGACACACA

4.5获取所有ebola病毒基因组的start序列

esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start 1 -seq_stop 10 > starts.fa

$ cat starts.fa |head
>KR105345.1:1-10 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G6089.1, partial genome
ATAATTTTCC
>KR105328.1:1-10 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G5844.1, partial genome
GATTAATAAT
>KR105323.1:1-10 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G5743.1, partial genome
GATTAATAAT
>KR105302.1:1-10 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G5295.1, partial genome
GATTAATAAT
>KR105295.1:1-10 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G5039.1, partial genome
ATTAATAATT
···

4.6可以把上面命令写入脚本

假如名字为get_seq.sh
内容如下

esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start 1 -seq_stop 10

则可以直接运行

bash get_seq.sh > starts.fa

5 查看quality和起始密码等具体信息

5.1看前 1 W行中的质量差的数据数目

$ cat SRR1553605_1.fastq |head -10000|grep '###'|wc -l
      82

5.2 看前1w行中的ATG

cat SRR1553605_1.fastq |head -10000|grep 'ATG' --color=always|head
ATG

5.3 看有无类似TATAbox

$ cat SRR1553605_1.fastq|head -10000|grep 'TATA' --color=always|head
TATA

5.4看有无类似illumina接头序列GATCGGA

 cat SRR1553605_1.fastq|head -10000|grep 'GATCGGA' --color=always|head

6 fastqc质量控制解释结果

chmod +x fastqc 
mkdir -p ~/bin
echo 'export PATH=~/bin:$PATH' >> ~/.profile
source ~/.profile
ln -s ~/src/FastQC/fastqc ~/bin/fastqc
fastqc -h
fastqc *.fastq
质量差
质量好

7碱基质量矫正base quality trimming

借用碱基矿工的这部分内容

当我们理解了fq数据之后，做这些过滤就不会很难，你也完全可以自己编写工具来进行个性化的过滤。目前也已有很多工具用来切除接头序列和低质量碱基，比如SOAPnuke、cutadapt、untrimmed等不下十个，但这其中比较方便好用的是Trimmomatic（也是一个java程序）、sickle和seqtk。Trimmomatic的好处在于，它不但可以用来切除illumina测序平台的接头序列，还可以去除由我们自己指定的特定接头序列，而且同时也能够过滤read末尾的低质量序列，sickle和seqtk只能去除低质量碱基。具体的原理就是通过滑动一定长度的窗口，计算窗口内的碱基平均质量，如果过低，就直接往后全部切除，注！意！不是挖掉read中的这部分低质量序列，而是像切菜一样，直接从低质量区域开始把这条read后面的所有其它碱基全！部！剁！掉！否则就是在人为改变实际的基因组序列情况。

7.1安装Trimmomatic

curl -O http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip
cd Trimmomatic-0.39/
java -jar trimmomatic-0.39.jar

看是否可用

Usage: 
       PE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] [-validatePairs] [-basein <inputBase> | <inputFile1> <inputFile2>] [-baseout <outputBase> | <outputFile1P> <outputFile1U> <outputFile2P> <outputFile2U>] <trimmer1>...
   or: 
       SE [-version] [-threads <threads>] [-phred33|-phred64] [-trimlog <trimLogFile>] [-summary <statsSummaryFile>] [-quiet] <inputFile> <outputFile> <trimmer1>...
   or: 
       -version

echo '#!/bin/bash' > ~/bin/trimmomatic
echo 'java -jar ~/src/Trimmomatic-0.39/trimmomatic-0.39.jar $@' >> ~/bin/trimmomatic
chmod +x ~/bin/trimmomatic
trimmomatic

7.2 PE命令如下

trimmomatic PE SRR1553610_1.fastq SRR1553610_2.fastq trim_1P.fq trim_1U.fq trim_2P.fq trim_2U.fq ILLUMINACLIP:TruSeq3-PE-2.fa:2:10:10 SLIDINGWINDOW:4:30 TRAILING:20

关于trimmomatic相关内容请见https://zhuanlan.zhihu.com/p/28924793
https://www.jianshu.com/p/a8935adebaae
结果会产生以下四个文件

trim_1P.fq
trim_1U.fq
trim_2P.fq
trim_2U.fq

8 模式匹配和adapter trimming

8.1 简单的正则表达式匹配

#ATG开头
  768  cat SRR1553605_1.fastq |egrep "^ATG" --color=always|head
#ATG结尾  
769  cat SRR1553605_1.fastq |egrep "ATG$" --color=always|head
#ATG结尾  
770  cat SRR1553605_1.fastq |egrep "ATG\$" --color=always|head
# TAAATA或TACCTA
 771  cat SRR1553605_1.fastq |egrep "TA(AA,CC)TA" --color=always|head
# TAAAA或TATAA 
 772  cat SRR1553605_1.fastq |egrep "TA[A,T]AA" --color=always|head
# TAAAA  TATAA
  773  cat SRR1553605_1.fastq |egrep "TA(A,T)AA" --color=always|head
#TA和AA中间有0或多个A
  774  cat SRR1553605_1.fastq |egrep "TA(A*)AA" --color=always|head
#TA和TA之间有0个或多个A
  775  cat SRR1553605_1.fastq |egrep "TA(A*)TA" --color=always|head
#TA和TA中间有1个或多个A
  776  cat SRR1553605_1.fastq |egrep "TA(A+)TA" --color=always|head
#TAA和TA中间有2到5个A
  777  cat SRR1553605_1.fastq |egrep "TAA{2,5}TA" --color=always|head

8.2 产生一个adapter文件

echo ">adapter" > adapter.fa
echo "AGATCGGAAGAG">>adapter.fa

8.3为trimmatic的adaptor文件夹创建软连接

ln -s ~/src/Trimmomatic-0.39/adapters/TruSeq3-PE.fa
ln -s ~/src/Trimmomatic-0.39/adapters/TruSeq3-PE-2.fa
ln -s ~/src/Trimmomatic-0.39/adapters/TruSeq3-SE.fa

9安装使用Blast

9.1 安装BLAST+

cd ~/src
curl -O ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.9.0+-x64-macosx.tar.gz
tar zxvf ncbi-blast-2.9.0+-x64-macosx.tar.gz 
echo 'export PATH=$PATH:~/src/ncbi-blast-2.9.0+/bin' >> ~/.profile
# for linux
#echo 'export PATH=$PATH:~/src/ncbi-blast-2.9.0+/bin' >> ~/.bashrc
source ~/.profile
blastn -h

已经可用

9.2 blast几个基本问题

1 我们想发现什么？query sequence
2 到哪里寻找？数据库database也就是target sequence
3 如何寻找？ search type

9.3 make一个blast 数据库

建一个Ebola病毒的基因组序列，因为index的时候会产生很多文件，所以建立一个新文件夹，命名为refs
因为reference可能包含很多fasta records，所以把它转变成blast database，以便程序可以处理

cd refs
makeblastdb -in KM233090.fa -dbtype nucl

会显示如下


Building a new DB, current time: 06/23/2019 19:10:50
New DB name:   /Users/ucco/refs/KM233090.fa
New DB title:  KM233090.fa
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 1 sequences in 0.000593901 seconds.
ucco:refs ucco$ ls
KM233090.fa KM233090.fa.nhr KM233090.fa.nin KM233090.fa.nsq KM233118.fa

然后会产生以下文件

-rw-r--r--  1 ucco  staff   162B  6 23 19:10 KM233090.fa.nhr
-rw-r--r--  1 ucco  staff    88B  6 23 19:10 KM233090.fa.nin
-rw-r--r--  1 ucco  staff   4.6K  6 23 19:10 KM233090.fa.nsq

9.4建立一个query序列

head -2 KM233090.fa >query.fa
cat query.fa 

>KM233090.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3816, complete genome
AAATTGTTACTGTAATCATACCTGGTTTGTTTCAGAGCCATATCACCAAGATAGAGAACAACCTAGGTCT

9.5 搜寻query序列（在nucleotide database）

9.5.1显示逐个碱基配对

blastn -db KM233090.fa -query query.fa|more

>KM233090.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3816, 
complete genome
Length=18799

 Score = 130 bits (70),  Expect = 6e-34
 Identities = 70/70 (100%), Gaps = 0/70 (0%)
 Strand=Plus/Plus

Query  1   AAATTGTTACTGTAATCATACCTGGTTTGTTTCAGAGCCATATCACCAAGATAGAGAACA  60
           ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  1   AAATTGTTACTGTAATCATACCTGGTTTGTTTCAGAGCCATATCACCAAGATAGAGAACA  60

Query  61  ACCTAGGTCT  70
           ||||||||||
Sbjct  61  ACCTAGGTCT  70

9.5.2 不需要显示base by base alignment

blastn -db KM233090.fa -query query.fa -outfmt 6

KM233090.1  KM233090.1  100.000 70  0   0   1   70  1   70

9.5.3 写入header information

blastn -db KM233090.fa -query query.fa -outfmt 7
# BLASTN 2.9.0+
# Query: KM233090.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-G3816, complete genome
# Database: KM233090.fa
# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 1 hits found
KM233090.1  KM233090.1  100.000 70  0   0   1   70  1   70  6.02e-34    130
# BLAST processed 1 queries

9.5.6其他形式

blastn -db KM233090.fa -query query.fa -outfmt '6 qseqid sseqid qlen length mismatch'
KM233090.1  KM233090.1  70  70  0

9.6 获取目标database的所有genome

esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta > all-genomes.fa

9.7 制作所有以上genome的blast database

makeblastdb -in all-genomes.fa -dbtype nucl

Building a new DB, current time: 06/23/2019 20:39:06
New DB name:   /Users/ucco/refs/all-genomes.fa
New DB title:  all-genomes.fa
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 249 sequences in 0.0657458 seconds.

9.8获取genome的特定区域序列

efetch -db nucleotide -id KM233118 -format fasta -seq_start 1 -seq_stop 1000 > 1000.fa

cat 1000.fa |head
>KM233118.1:1-1000 Zaire ebolavirus isolate Ebola virus/H.sapiens-wt/SLE/2014/Makona-NM042.3, complete genome
AATCATACCTGGTTTGTTTCAGAGCCATATCACCAAGATAGAGAACAACCTAGGTCTCCGGAGGGGGCAA
GGGCATCAGTGTGCTCAGTTGAAAATCCCTTGTCAACATCTAGGCCTTATCACATCACAAGTTCCGCCTT
AAACTCTGCAGGGTGATCCAACAACCTTAATAGCAACATTATTGTTAAAGGACAGCATTAGTTCACAGTC
AAACAAGCAAGATTGAGAATTAACTTTGATTTTGAACCTGAACACCCAGAGGACTGGAGACTCAACAACC
CTAAAGCCTGGGGTAAAACATTAGAAATAGTTTAAAGACAAATTGCTCGGAATCACAAAATTCCGAGTAT
GGATTCTCGTCCTCAGAAAGTCTGGATGACGCCGAGTCTCACTGAATCTGACATGGATTACCACAAGATC
TTGACAGCAGGTCTGTCCGTTCAACAGGGGATTGTTCGGCAAAGAGTCATCCCAGTGTATCAAGTAAACA
ATCTTGAGGAAATTTGCCAACTTATCATACAGGCCTTTGAAGCTGGTGTTGATTTTCAAGAGAGTGCGGA
CAGTTTCCTTCTCATGCTTTGTCTTCATCATGCGTACCAAGGAGATTACAAACTTTTCTTGGAAAGTGGC

10 blastdbcmd, blastp, blastx, tblastn的用法

1 将新病毒和1999年爆发的病毒进行比较
https://www.ncbi.nlm.nih.gov/genome/genomes/4887?
2 RefSeq accession number NC_002549.1, Nucleoprotein example: NP_066243.1

10.1获取feature table

关于feature table请看上传基因组数据到NCBI中的详细解释

efetch -db nucleotide -id NC_002549  --format fasta > refs/NC_002549.fa
esearch -db protein -query PRJNA14703 | efetch -format ft

注意，这里PRJNA14703已经不存在了，应该是重新进行了提交申请。



所以我下载了PRJNA485481，比较大200多M

fetch -db nucleotide -id NC_002549  --format fasta > refs/NC_002549.fa
esearch -db protein -query PRJNA485481 | efetch -format fasta > NC_002549-prot.fa

10.2 为核苷酸和蛋白质制作Blast 数据

makeblastdb -in refs/NC_002549.fa -dbtype nucl -parse_seqids
makeblastdb -in refs/NC_002549-prot.fa -dbtype prot -parse_seqids

10.3 blastdbcmd的多种用法

blastdbcmd有很多参数可以查询和格式化blast database的内容
列出database中的所有内容

blastdbcmd -db refs/NC_002549.fa -entry 'all' | head -3

>NC_002549.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-tc/COD/1976/Yambuku-Mayinga, complete genome
CGGACACACAAAAAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATAACTATGAGGAAGATTAATAATTTTCCTCTC
ATTGAAATTTATATCGGAATTTAAATTGAAATTGTTACTGTAATCACACCTGGTTTGTTTCAGAGCCACATCACAAAGAT

获取database中一系列的核苷酸entry

blastdbcmd -db refs/NC_002549.fa -entry 'all' -outfmt '%a %s' -range 1-10 -strand minus

NC_002549.1 TGTGTGTCCG

输出特定格式的database内容

blastdbcmd -db refs/NC_002549.fa -entry 'all' -outfmt '%a %l'

NC_002549.1 18959

其中一些输出格式的参数表示如下

        %f means sequence in FASTA format
        %s means sequence data (without defline)
        %a means accession
        %g means gi
        %o means ordinal id (OID)
        %i means sequence id
        %t means sequence title
        %l means sequence length
        %h means sequence hash value
        %T means taxid
        %X means leaf-node taxids
        %e means membership integer
        %L means common taxonomic name
        %C means common taxonomic names for leaf-node taxids
        %S means scientific name
        %N means scientific names for leaf-node taxids
        %B means BLAST name
        %K means taxonomic super kingdom
        %P means PIG
        %m means sequence masking data

列出每个蛋白的长度

blastdbcmd -db refs/NC_002549-prot.fa -entry 'all' -outfmt '%a %l'|head
YP_009646896.1 315
YP_009646895.1 64
YP_009646894.1 138
YP_009646893.1 326
YP_009646892.1 172
YP_009646891.1 227
YP_009646890.1 55
YP_009646889.1 68
YP_009646888.1 91
YP_009646887.1 90

提取一个蛋白entry到文件

blastdbcmd -db refs/NC_002549-prot.fa -entry 'NP_066243.1' > query-p.fa
cat query-p.fa

>NP_066243.1 nucleoprotein [Zaire ebolavirus]
MDSRPQKIWMAPSLTESDMDYHKILTAGLSVQQGIVRQRVIPVYQVNNLEEICQLIIQAFEAGVDFQESADSFLLMLCLH
HAYQGDYKLFLESGAVKYLEGHGFRFEVKKRDGVKRLEELLPAVSSGKNIKRTLAAMPEEETTEANAGQFLSFASLFLPK
LVVGEKACLEKVQRQIQVHAEQGLIQYPTAWQSVGHMMVIFRLMRTNFLIKFLLIHQGMHMVAGHDANDAVISNSVAQAR
FSGLLIVKTVLDHILQKTERGVRLHPLARTAKVKNEVNSFKAALSSLAKHGEYAPFARLLNLSGVNNLEHGLFPQLSAIA
LGVATAHGSTLAGVNVGEQYQQLREAATEAEKQLQQYAESRELDHLGLDDQEKKILMNFHQKKNEISFQQTNAMVTLRKE
RLAKLTEAITAASLPKTSGHYDDDDDIPFPGPINDDDNPGHQDDDPTDSQDTTIPDVVVDPDDGSYGEYQSYSENGMNAP
DDLVLFDLDEDDEDTKPVPNRSTKGGQQKNSQKGQHIEGRQTQSRPIQNVPGPHRTIHHASAPLTDNDRRNEPSGSTSPR
MLTPINEEADPLDDADDETSSLPPLESDDEEQDRDGTSNRTPTVAPPAPVYRDHSEKKELPQDEQQDQDHTQEARNQDSD
NTQSEHSFEEMYRHILRSQGPFDAVLYYHMMKDEPVVFSTSDGKEYTYPDSLEEEYPPWLTEKEAMNEENRFVTLDGQQF
YWPVMNHKNKFMAILQHHQ

10.4 运行blastp

blastp -db refs/NC_002549-prot.fa -query query-p.fa

10.5blastp

先在genebank获取一个蛋白编码序列的区域

blastdbcmd -db refs/NC_002549.fa -entry 'NC_002549' -range 470-2689  > nucleotide.fa
cat nucleotide.fa |head

>NC_002549.1:470-2689 Zaire ebolavirus isolate Ebola virus/H.sapiens-tc/COD/1976/Yambuku-Mayinga, complete genome
ATGGATTCTCGTCCTCAGAAAATCTGGATGGCGCCGAGTCTCACTGAATCTGACATGGATTACCACAAGATCTTGACAGC
AGGTCTGTCCGTTCAACAGGGGATTGTTCGGCAAAGAGTCATCCCAGTGTATCAAGTAAACAATCTTGAAGAAATTTGCC
AACTTATCATACAGGCCTTTGAAGCAGGTGTTGATTTTCAAGAGAGTGCGGACAGTTTCCTTCTCATGCTTTGTCTTCAT
CATGCGTACCAGGGAGATTACAAACTTTTCTTGGAAAGTGGCGCAGTCAAGTATTTGGAAGGGCACGGGTTCCGTTTTGA
AGTCAAGAAGCGTGATGGAGTGAAGCGCCTTGAGGAATTGCTGCCAGCAGTATCTAGTGGAAAAAACATTAAGAGAACAC
TTGCTGCCATGCCGGAAGAGGAGACAACTGAAGCTAATGCCGGTCAGTTTCTCTCCTTTGCAAGTCTATTCCTTCCGAAA
TTGGTAGTAGGAGAAAAGGCTTGCCTTGAGAAGGTTCAAAGGCAAATTCAAGTACATGCAGAGCAAGGACTGATACAATA
TCCAACAGCTTGGCAATCAGTAGGACACATGATGGTGATTTTCCGTTTGATGCGAACAAATTTTCTGATCAAATTTCTCC
TAATACACCAAGGGATGCACATGGTTGCCGGGCATGATGCCAACGATGCTGTGATTTCAAATTCAGTGGCTCAAGCTCGT

把核酸序列比对蛋白database

blastx -db refs/NC_002549-prot.fa -query nucleotide.fa | more

部分结果如下

Sequences producing significant alignments:                          (Bits)     Value

NP_066243.1 nucleoprotein [Zaire ebolavirus]                          1327       0.0   
YP_003815432.1 nucleoprotein [Bundibugyo ebolavirus]                  1007       0.0   
YP_003815423.1 nucleoprotein [Tai Forest ebolavirus]                  785        0.0   
YP_009513274.1 nucleoprotein [Bombali ebolavirus]                     769        0.0   
NP_690580.1 nucleoprotein [Reston ebolavirus]                         769        0.0   
YP_138520.1 nucleoprotein [Sudan ebolavirus]                          750        0.0   
YP_004928135.1 nucleoprotein [Lloviu cuevavirus]                      576        0.0

10.5 获取核蛋白序列

efetch -db protein -id NP_066243.1 -format fasta > NP_066243.fa

cat NP_066243.fa 
>NP_066243.1 nucleoprotein [Zaire ebolavirus]
MDSRPQKIWMAPSLTESDMDYHKILTAGLSVQQGIVRQRVIPVYQVNNLEEICQLIIQAFEAGVDFQESA
DSFLLMLCLHHAYQGDYKLFLESGAVKYLEGHGFRFEVKKRDGVKRLEELLPAVSSGKNIKRTLAAMPEE
ETTEANAGQFLSFASLFLPKLVVGEKACLEKVQRQIQVHAEQGLIQYPTAWQSVGHMMVIFRLMRTNFLI
KFLLIHQGMHMVAGHDANDAVISNSVAQARFSGLLIVKTVLDHILQKTERGVRLHPLARTAKVKNEVNSF
KAALSSLAKHGEYAPFARLLNLSGVNNLEHGLFPQLSAIALGVATAHGSTLAGVNVGEQYQQLREAATEA
EKQLQQYAESRELDHLGLDDQEKKILMNFHQKKNEISFQQTNAMVTLRKERLAKLTEAITAASLPKTSGH
YDDDDDIPFPGPINDDDNPGHQDDDPTDSQDTTIPDVVVDPDDGSYGEYQSYSENGMNAPDDLVLFDLDE
DDEDTKPVPNRSTKGGQQKNSQKGQHIEGRQTQSRPIQNVPGPHRTIHHASAPLTDNDRRNEPSGSTSPR
MLTPINEEADPLDDADDETSSLPPLESDDEEQDRDGTSNRTPTVAPPAPVYRDHSEKKELPQDEQQDQDH
TQEARNQDSDNTQSEHSFEEMYRHILRSQGPFDAVLYYHMMKDEPVVFSTSDGKEYTYPDSLEEEYPPWL
TEKEAMNEENRFVTLDGQQFYWPVMNHKNKFMAILQHHQ

tblastn -db refs/NC_002549.fa -query NP_066243.fa  | more

部分结果如下


Query= NP_066243.1 nucleoprotein [Zaire ebolavirus]

Length=739
                                                                      Score        E
Sequences producing significant alignments:                          (Bits)     Value

NC_002549.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-tc/COD...  1328       0.0  


>NC_002549.1 Zaire ebolavirus isolate Ebola virus/H.sapiens-tc/COD/1976/Yambuku-Mayinga, 
complete genome
Length=18959

 Score = 1328 bits (3436),  Expect = 0.0, Method: Compositional matrix adjust.
 Identities = 739/739 (100%), Positives = 739/739 (100%), Gaps = 0/739 (0%)
 Frame = +2

Query  1     MDSRPQKIWMAPSLTESDMDYHKILTAGLSVQQGIVRQRVIPVYQVNNLEEICQLIIQAF  60
             MDSRPQKIWMAPSLTESDMDYHKILTAGLSVQQGIVRQRVIPVYQVNNLEEICQLIIQAF
Sbjct  470   MDSRPQKIWMAPSLTESDMDYHKILTAGLSVQQGIVRQRVIPVYQVNNLEEICQLIIQAF  649

10.6 获取2014 ebola project的所有蛋白数据

esearch -db protein -query PRJNA257197 | efetch -format fasta > refs/ebola-2014.fa
makeblastdb -in refs/ebola-2014.fa -dbtye prot -parse_seqids

cat NP_066243.fa 
>NP_066243.1 nucleoprotein [Zaire ebolavirus]
MDSRPQKIWMAPSLTESDMDYHKILTAGLSVQQGIVRQRVIPVYQVNNLEEICQLIIQAFEAGVDFQESA
DSFLLMLCLHHAYQGDYKLFLESGAVKYLEGHGFRFEVKKRDGVKRLEELLPAVSSGKNIKRTLAAMPEE
ETTEANAGQFLSFASLFLPKLVVGEKACLEKVQRQIQVHAEQGLIQYPTAWQSVGHMMVIFRLMRTNFLI
KFLLIHQGMHMVAGHDANDAVISNSVAQARFSGLLIVKTVLDHILQKTERGVRLHPLARTAKVKNEVNSF
KAALSSLAKHGEYAPFARLLNLSGVNNLEHGLFPQLSAIALGVATAHGSTLAGVNVGEQYQQLREAATEA
EKQLQQYAESRELDHLGLDDQEKKILMNFHQKKNEISFQQTNAMVTLRKERLAKLTEAITAASLPKTSGH
YDDDDDIPFPGPINDDDNPGHQDDDPTDSQDTTIPDVVVDPDDGSYGEYQSYSENGMNAPDDLVLFDLDE
DDEDTKPVPNRSTKGGQQKNSQKGQHIEGRQTQSRPIQNVPGPHRTIHHASAPLTDNDRRNEPSGSTSPR
MLTPINEEADPLDDADDETSSLPPLESDDEEQDRDGTSNRTPTVAPPAPVYRDHSEKKELPQDEQQDQDH
TQEARNQDSDNTQSEHSFEEMYRHILRSQGPFDAVLYYHMMKDEPVVFSTSDGKEYTYPDSLEEEYPPWL
TEKEAMNEENRFVTLDGQQFYWPVMNHKNKFMAILQHHQ
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 217,826评论 6赞 506
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 92,968评论 3赞 395
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 164,234评论 0赞 354
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 58,562评论 1赞 293
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 67,611评论 6赞 392
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 51,482评论 1赞 302
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 40,271评论 3赞 418
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 39,166评论 0赞 276
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 45,608评论 1赞 314
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,814评论 3赞 336
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,926评论 1赞 348
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 35,644评论 5赞 346
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 41,249评论 3赞 329
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,866评论 0赞 22
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,991评论 1赞 269
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 48,063评论 3赞 370
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 44,871评论 2赞 354

推荐阅读更多精彩内容