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

Y大宽

发布于 2019-06-24 10:48:39

2.3K0

发布于 2019-06-24 10:48:39

文章被收录于专栏：Y大宽Y大宽

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

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 删除

https

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

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

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度

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

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

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

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐