通过简单数据熟悉Linux下生物信息学各种操作
2019-06-23 本文已影响38人
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
其他几种下载方式，看
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
