前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >通过简单数据熟悉Linux下生物信息学各种操作

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

作者头像
Y大宽
发布2019-06-24 10:48:39
2.4K0
发布2019-06-24 10:48:39
举报
文章被收录于专栏:Y大宽

原地址

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

代码语言:javascript
复制
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

代码语言:javascript
复制
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下载数据

代码语言:javascript
复制
prefetch SRR1553610
find ~/ncbi

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

代码语言:javascript
复制
/Users/ucco/ncbi
/Users/ucco/ncbi/public
/Users/ucco/ncbi/public/sra
/Users/ucco/ncbi/public/sra/SRR1553610.sra

sra格式转变为fq格式

代码语言:javascript
复制
fastq-dump --split-files SRR1553610.sra
代码语言:javascript
复制
wc -l *.fastq
  879348 SRR1553610_1.fastq
  879348 SRR1553610_2.fastq
 1758696 total

一次下载多个文件

代码语言:javascript
复制
$ echo SRR1553608 > sra.ids
$ echo SRR1553605 >> sra.ids
$ prefetch --option-file sra.ids

其他几种下载方式,看

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

4 通过EDirect获取序列

4.1根据locus获取序列

代码语言:javascript
复制
efetch -db nucleotide -id KM233090 -format fasta > KM233090.fa

4.2 根据accession number获取序列

代码语言:javascript
复制
efetch -db nucleotide -id 667853062 -format fasta > 667853062.fa
代码语言:javascript
复制
ucco:yeast ucco$ cat KM233090.fa |wc
     270     277   19169
ucco:yeast ucco$ cat 667853062.fa |wc
     270     277   19169

4.3同时获取多个序列

代码语言:javascript
复制
efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta > multi.fa

看下header files

代码语言:javascript
复制
$ 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只获取其中一部分

代码语言:javascript
复制
efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta -seq_start 1 -seq_stop 10 > multi.fa
代码语言:javascript
复制
$ 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序列

代码语言:javascript
复制
esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start 1 -seq_stop 10 > starts.fa
代码语言:javascript
复制
$ 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

内容如下

代码语言:javascript
复制
esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start 1 -seq_stop 10

则可以直接运行

代码语言:javascript
复制
bash get_seq.sh > starts.fa

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

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

代码语言:javascript
复制
$ cat SRR1553605_1.fastq |head -10000|grep '###'|wc -l
      82

5.2 看前1w行中的ATG

代码语言:javascript
复制
cat SRR1553605_1.fastq |head -10000|grep 'ATG' --color=always|head

ATG

5.3 看有无类似TATAbox

代码语言:javascript
复制
$ cat SRR1553605_1.fastq|head -10000|grep 'TATA' --color=always|head

TATA

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

代码语言:javascript
复制
 cat SRR1553605_1.fastq|head -10000|grep 'GATCGGA' --color=always|head

6 fastqc质量控制解释结果

代码语言:javascript
复制
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

代码语言:javascript
复制
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

看是否可用

代码语言:javascript
复制
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
代码语言:javascript
复制
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命令如下

代码语言:javascript
复制
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

结果会产生以下四个文件

代码语言:javascript
复制
trim_1P.fq
trim_1U.fq
trim_2P.fq
trim_2U.fq

8 模式匹配和adapter trimming

8.1 简单的正则表达式匹配

代码语言:javascript
复制
#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文件

代码语言:javascript
复制
echo ">adapter" > adapter.fa
echo "AGATCGGAAGAG">>adapter.fa

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

代码语言:javascript
复制
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+

代码语言:javascript
复制
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,以便程序可以处理

代码语言:javascript
复制
cd refs
makeblastdb -in KM233090.fa -dbtype nucl

会显示如下

代码语言:javascript
复制
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

然后会产生以下文件

代码语言:javascript
复制
-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序列

代码语言:javascript
复制
head -2 KM233090.fa >query.fa
cat query.fa 
代码语言:javascript
复制
>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显示逐个碱基配对

代码语言:javascript
复制
blastn -db KM233090.fa -query query.fa|more
代码语言:javascript
复制
>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

代码语言:javascript
复制
blastn -db KM233090.fa -query query.fa -outfmt 6
代码语言:javascript
复制
KM233090.1  KM233090.1  100.000 70  0   0   1   70  1   70

9.5.3 写入header information

代码语言:javascript
复制
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其他形式

代码语言:javascript
复制
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

代码语言:javascript
复制
esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta > all-genomes.fa

9.7 制作所有以上genome的blast database

代码语言:javascript
复制
makeblastdb -in all-genomes.fa -dbtype nucl
代码语言:javascript
复制
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的特定区域序列

代码语言:javascript
复制
efetch -db nucleotide -id KM233118 -format fasta -seq_start 1 -seq_stop 1000 > 1000.fa
代码语言:javascript
复制
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

代码语言:javascript
复制
efetch -db nucleotide -id NC_002549  --format fasta > refs/NC_002549.fa
esearch -db protein -query PRJNA14703 | efetch -format ft
本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
原始发表:2019.06.23 ,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1下载酵母基因组gff格式文件
  • 2安装sratoolkit
  • 3下载数据
  • 其他几种下载方式,看
  • 4 通过EDirect获取序列
    • 4.1根据locus获取序列
      • 4.2 根据accession number获取序列
        • 4.3同时获取多个序列
          • 4.4只获取其中一部分
          • 4.5获取所有ebola病毒基因组的start序列
          • 4.6可以把上面命令写入脚本
          • 5 查看quality和起始密码等具体信息
            • 5.1看前 1 W行中的质量差的数据数目
              • 5.2 看前1w行中的ATG
                • 5.3 看有无类似TATAbox
                  • 5.4看有无类似illumina接头序列GATCGGA
                  • 6 fastqc质量控制解释结果
                  • 7碱基质量矫正base quality trimming
                    • 7.1安装Trimmomatic
                      • 7.2 PE命令如下
                      • 8 模式匹配和adapter trimming
                        • 8.1 简单的正则表达式匹配
                          • 8.2 产生一个adapter文件
                            • 8.3为trimmatic的adaptor文件夹创建软连接
                            • 9安装使用Blast
                              • 9.1 安装BLAST+
                                • 9.2 blast几个基本问题
                                  • 9.3 make一个blast 数据库
                                    • 9.4建立一个query序列
                                      • 9.5 搜寻query序列(在nucleotide database)
                                        • 9.5.1显示逐个碱基配对
                                        • 9.5.2 不需要显示base by base alignment
                                        • 9.5.3 写入header information
                                        • 9.5.6其他形式
                                      • 9.6 获取目标database的所有genome
                                        • 9.7 制作所有以上genome的blast database
                                          • 9.8获取genome的特定区域序列
                                          • 10 blastdbcmd, blastp, blastx, tblastn的用法
                                            • 10.1获取feature table
                                            相关产品与服务
                                            数据库
                                            云数据库为企业提供了完善的关系型数据库、非关系型数据库、分析型数据库和数据库生态工具。您可以通过产品选择和组合搭建,轻松实现高可靠、高可用性、高性能等数据库需求。云数据库服务也可大幅减少您的运维工作量,更专注于业务发展,让企业一站式享受数据上云及分布式架构的技术红利!
                                            领券
                                            问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档