SeqKit的学习 --20191017
软件的介绍
SeqKit是一种跨平台的、极快的,全面的fasta/q处理工具。SeqKit为所有的主流操作系统提供了一种可执行的双元文件,包括Windows,Linux,Mac OS X,并且不依赖于任何的配置或预先配置就可以直接使用。
软件的安装
## Install via conda
conda install -c bioconda seqkit
软件的命令
## 序列和子序列
**seq** 转换序列(序列颠倒,序列互补,提取ID)
**subseq** 从区域/gtf/bed中获得序列,包括侧面的序列
**sliding** 滑动序列,支持环式基因组
**stats** 对FASTA/Q files进行简单统计
**faidx** 创造fasta索引文件并提取子序列
**watch** 检测并连线序列特点的柱状图
**sana** 清除质量不好的单线的fastq文件
## 格式转换
**fx2tab** 将FASTA/Q 文件转变成表格形式 (1th: name/ID, 2nd: sequence, 3rd: quality)
**tab2fx** 转变表格形式为fasta/q格式
**fq2fa** 转变fastq文件为fasta文件
**convert** 在Sanger, Solexa and Illumina中转换fastq的质量编码
**translate** 将DNA/RNA序列转变成蛋白序列(支持模棱两可的碱基)
## 搜索
**grep** 根据ID/名称/序列/序列motif 搜索序列,且允许错配
**locate** 定位子序列/motif,且允许错配
**fish** 使用本地比对在较大序列中寻找短序列
**amplicon** 经由引物检索扩增子(或它附近特定的区域)
## bam文件的处理和监视
**bam** 监视和连线bam文件记录特点的直方图
## 设置参数
**head** 打印第一个Nfasta/q的记录
**range** 在一个范围内(start:end)打印fasta/q的记录
**sample** 通过数量或比例来体验序列
**rmdup** 通过id/名称/序列 来去除复制的序列
**duplicate** 复制N次的序列
**common** 通过id/名称/序列 发现多条序列中共有的序列
**split** 通过id/seq region/size/parts (mainly for FASTA) 将序列劈开成文件
**split2** 将序列通过大小或部分 劈开成文件
## 编辑
**replace** 通过规律表达来代替名字或序列
**rename** 重新命名复制的ID
**restart** 为环状基因组重新设置起始位置
**concat** 从多个文件中经由相同的ID来连接序列
**mutate** 编辑序列(点突,插入,删除)
## 排序
**shuffle** 变换序列位置
**sort** 将序列经由id/name/sequence 进行排序
软件命令详解
Sequence ID
大部分的软件,包括seqkit默认将主导的非空格字母作为ID。
FASTA header | ID |
---|---|
>123456 gene name | 123456 |
>longname | longname |
>gi|110645304|ref|NC_002516.2| Pseudomona | gi|110645304|ref|NC_002516.2| |
举例说明软件如何使用
##下载参考序列,一个fastq文件,两个fasta文件
wget http://data.biostarhandbook.com/reads/duplicated-reads.fq.gz
wget ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.1.genomic.fna.gz
wget ftp://ftp.ncbi.nih.gov/refseq/release/viral/viral.1.protein.faa.gz
对fastq文件进行一个概括浏览
$ seqkit stat *.gz
file format type num_seqs sum_len min_len avg_len max_len
duplicated-reads.fq.gz FASTQ DNA 15,000 1,515,000 101 101 101
viral.1.1.genomic.fna.gz FASTA DNA 6 195,842 5,386 32,640.3 154,675
viral.1.protein.faa.gz FASTA Protein 112 55,855 38 498.7 3,122
在fasta/q文件中获取每条序列的GC含量
GC content
$ seqkit fx2tab --name --only-id --gc viral*.fna.gz
NC_001798.2 70.38
NC_030692.1 50.10
NC_027892.1 40.57
NC_029642.1 39.88
NC_001607.1 50.07
NC_001422.1 44.76
Custom bases (A, C and A+C)
$ seqkit fx2tab -H -n -i -B a -B c -B ac viral.1.1.genomic.fna.gz
#name seq qual a c ac
NC_001798.2 14.87 35.03 49.90
NC_030692.1 25.03 25.25 50.28
NC_027892.1 29.68 19.44 49.11
NC_029642.1 30.14 19.22 49.36
NC_001607.1 25.30 25.54 50.84
NC_001422.1 23.97 21.48 45.45
附seqkit fx2tab 的使用(用来统计序列的信息)
Usage: seqkit fx2tab [flags]
Flags:
-a, --alphabet 打印字母表字母
-q, --avg-qual 打印read的平均质量
-B, --base-content strings 输出指定的碱基含量
-g, --gc 输出GC含量
-H, --header-line 打印标题列
-l, --length 打印序列的长度
-n, --name 只打印名字,而没有序列或者质量
-i, --only-id 只打印基因的ID
从fastq/a中根据名字和ID提取序列子集
$ seqkit sample --proportion 0.001 duplicated-reads.fq.gz \
| seqkit seq --name --only-id > id.txt ##管道符前面的命令是随机取总文件0.1%的序列,管道符后面的是提取前面的符合要求的序列的ID。
## 查看ID list内容
$ head id.txt
SRR1972739.2996
SRR1972739.3044
SRR1972739.3562
##通过ID list文件来搜索序列
$ seqkit grep --pattern-file id.txt duplicated-reads.fq.gz > duplicated-reads.subset.fq.gz
seqkit sample 的使用
-n, --number int 通过数量随机提取序列(结果也许并不完全匹配)
-p, --proportion float 通过比例随机提取序列
-s, --rand-seed int 随机函数
-2, --two-pass 减少内存占用
举例:随机抽取序列
seqkit sample -n 10000 -s 11 test1_1.fq -o sample.fq
seqkit sample -p 0.1 -s 11 test1_1.fq -o sample.fq
从fasta/q序列中找到合并碱基并找到它的位置(这个仿佛有点难度,不报错也不打印内容到屏幕)??
$ seqkit fx2tab -n -i -a viral.1.1.genomic.fna.gz \
| csvtk -H -t grep -f 4 -r -i -p "[^ACGT]"
## 定位简并碱基,e.g, N and K,这个仿佛也有问题,不报错也不打印内容到屏幕???
seqkit grep --pattern-file id2.txt viral.1.1.genomic.fna.gz \
| seqkit locate --ignore-case --only-positive-strand --pattern K+ --pattern N+
移去相同序列中重复的fasta/q记录
$ seqkit rmdup --by-seq --ignore-case duplicated-reads.fq.gz > duplicated-reads.uniq.fq.gz
[INFO] 7172 duplicated records removed
seqkit rmdup 的使用(用来通过id/名称/序列来去除重复的序列)
Usage: seqkit rmdup [flags]
Flags:
-n, --by-name 通过全名而不是id
-s, --by-seq 通过序列
-D, --dup-num-file string 保存数量的文件并列出重复的seqs
-d, --dup-seqs-file string 保存重复seqs的文件
-i, --ignore-case 忽视字母大小写
在fastq/a序列中定位motif/子序列/酶切位点
$ cat enzymes.fa
>EcoRI
GAATTC
>MmeI
TCCRAC
>SacI
GAGCTC
>XcmI
CCANNNNNNNNNTGG
$ seqkit locate --degenerate --ignore-case --pattern-file enzymes.fa viral.1.1.genomic.fna.gz
输出的内容
seqID patternName pattern strand start end matched
gi|526245010|ref|NC_021865.1| MmeI TCCRAC + 1816 1821 TCCGAC
gi|526245010|ref|NC_021865.1| SacI GAGCTC + 19506 19511 GAGCTC
gi|526245010|ref|NC_021865.1| XcmI CCANNNNNNNNNTGG + 2221 2235 CCATATTTAGTGTGG
seqkit locate 的使用
Usage:
seqkit locate [flags]
-d, --degenerate 包含简并碱基模式和motif
--gtf 输出为GTF格式
-i, --ignore-case 忽视字母大小写
-m, --max-mismatch int 通过序列匹配时允许的最大错配
-G, --non-greedy 非贪婪模式,更快但是可能错过与其他重叠的motif
-P, --only-positive-strand 只搜索正链
-f, --pattern-file 模式或motif文件(fasta格式)
-p, --pattern strings 搜索pattern或motif
seqkit locate -i -d -p AUGGACUN test.fa
输出结果:
seqID patternName pattern strand start end matched
cel-mir-58a AUGGACUN AUGGACUN + 81 88 AUGGACUG
ath-MIR163 AUGGACUN AUGGACUN - 122 129 AUGGACUC
怎样通过长度来对大量的fasta文件进行排序
$ seqkit sort --by-length viral.1.1.genomic.fna.gz > viral.1.1.genomic.sorted.fa
seqkit sort 的使用
通过id/名称/序列/长度来排序,对于fasta格式的文件,使用flag -2 来减少内存的使用,不支持fastq格式的文件
Usage:seqkit sort [flags]
Flags:
-l, --by-length 通过序列长度
-n, --by-name 通过全名而不是id
-s, --by-seq 通过序列
-i, --ignore-case 忽视大小写
-r, --reverse 反向排序
-2, --two-pass 双流程模式读取文件来降低内存使用
根据标题信息来拆分fasta序列
## 先对fasta文件的标题进行概括浏览
$ seqkit head -n 3 viral.1.protein.faa.gz | seqkit seq --name --only-id
YP_009137150.1
YP_009137151.1
YP_009137152.1
## 将默认的ID转换成新的ID
$ seqkit head -n 3 viral.1.protein.faa.gz \
| seqkit seq --name --only-id --id-regexp "\[(.+)\]" 这一步并不懂???
Human alphaherpesvirus 2
Human alphaherpesvirus 2
Human alphaherpesvirus 2
## 根据header进行拆分,会生成一个文件夹
$ seqkit split --by-id --id-regexp "\[(.+)\]" viral.1.protein.faa.gz
seqkit head 的使用(首先打印第N条fasta/q记录)
Usage:seqkit head [flags]
Flags:
-n, --number int 先打印N个fasta/q的记录(默认是10)
seqkit seq 的使用
对序列进行转换(颠倒,互补,提取ID等)
Usage:seqkit seq [flags]
Flags:
-p, --complement 取互补序列
--dna2rna DNA到RNA
--rna2dna RNA到DNA
-G, --gap-letters string gap letters (default "- \t.")
-l, --lower-case 用小写字母打印序列
-M, --max-len int 只打印短于最大长度的的序列
-n, --name 只打印name
-g, --remove-gaps 移去组装序列中的gap
-r, --reverse 取反向序列
-i, --only-id 只打印ID而不是全名
-q, --qual 只打印品质???
-s, --seq 只打印序列
--id-regexp string 对于解析ID的正则表达式,(default "^(\\S+)\\s?")
seqkit split 的使用
usage:seqkit split [flags]
-i, --by-id 根据序列ID切割序列
-p, --by-part int 将一份文件分割成N份
-s, --by-size int 将一个文件按照N条序列进行分割
-O, --out-dir string 输出路径
-2, --two-pass 降低内存的使用
从一个text文件已知的字符串中搜索并替换fasta标题
$ seqkit grep --by-name --use-regexp --ignore-case --pattern hypothetical viral.1.protein.faa.gz \
| seqkit head -n 3 | seqkit seq --name ## 此命令也是有报错,说明没有执行好???
$ seqkit replace --kv-file changes.tsv --pattern "^([^ ]+ )(\w+) " \
--replacement "\${1}{kv} " --key-capt-idx 2 --keep-key viral.1.protein.faa.gz > renamed.fa
seqkit grep 的使用(通过ID/名称/序列/序列motif来搜索序列,允许错配)
Usage:seqkit grep [flags]
Examples:
1-based index 1 2 3 4 5 6 7 8 9 10
negative index 0-9-8-7-6-5-4-3-2-1
seq A C G T N a c g t n
1:1 A
2:4 C G T
-4:-2 c g t
-4:-1 c g t n
-1:-1 n
2:-2 C G T N a c g t
1:-1 A C G T N a c g t n
1:12 A C G T N a c g t n
-12:-1 A C G T N a c g t n
-n, --by-name 通过全名来匹配而不是ID
-s, --by-seq 搜索seq中的子seq,只搜索正链,通过-m/--max-mismatch 来允许错配
-d, --degenerate 包含简并碱基的pattern和motif(简并碱基:一个符号代替某两个或者更多碱基,如编译丙氨酸的可以有4个密码子:GCU\GCC\GCA\GCG,这时生物学上为了方便,用字母N指代UCAG四个碱基,故说编译丙氨酸的密码子是GCN,其中N就是简并碱基。)
-i, --ignore-case 忽视大小写
-v, --invert-match 输出不匹配此模式的内容
-m, --max-mismatch int 通过序列匹配时允许的最大错配
-p, --pattern strings 匹配模式,支持连续写多个模式,匹配任一模式即输出。
-f, --pattern-file string 支持匹配模式写到一个文件中,如要提取的序列ID。
-R, --region string 搜索特定区域序列,e.g 1:12 for first 12 bases, -12:-1 for last 12 bases
-r, --use-regexp 使用正则表达式,必须加入此参数,如^匹配首端。同-p联合使用。
举例:
seqkit grep -s -r -i -p ^atg cds.fa#选取正链中有起始密码子的序列
seqkit grep -f ID.txt test.fa > new.fa#根据ID提取序列
seqkit grep -s -d -i -p TTSAA#简并碱基使用。S 代表C or G.
seqkit grep -s -R 1:30 -i -r -p GCTGG##匹配限定到某区域
usage:seqkit replace(通过正则表达式来代替名字或序列)
-s, --by-seq 代替seq
-i, --ignore-case 忽视大小写
-K, --keep-key
-p, --pattern string 搜索正则表达式
-r, --replacement string 替换物
从两个配对端读数的文件提取配对的reads
## 首先创造两个不平衡的PE reads文件,整个过程不报错,但是不懂具体含义???
$ seqkit rmdup duplicated-reads.fq.gz \
| seqkit replace --pattern " .+" --replacement " 1" \
| seqkit sample --proportion 0.9 --rand-seed 1 --out-file read_1.fq.gz
$ seqkit rmdup duplicated-reads.fq.gz \
| seqkit replace --pattern " .+" --replacement " 2" \
| seqkit sample --proportion 0.9 --rand-seed 2 --out-file read_2.fq.gz
## overview
$ seqkit stat read_1.fq.gz read_2.fq.gz
file format type num_seqs sum_len min_len avg_len max_len
read_1.fq.gz FASTQ DNA 9,033 912,333 101 101 101
read_2.fq.gz FASTQ DNA 8,965 905,465 101 101 101
## sequence headers
$ seqkit head -n 3 read_1.fq.gz | seqkit seq --name
SRR1972739.1 1
SRR1972739.3 1
SRR1972739.4 1
$ seqkit head -n 3 read_2.fq.gz | seqkit seq --name
SRR1972739.1 2
SRR1972739.2 2
SRR1972739.3 2
## 首先提取两个文件序列的ID,并计算它们的交集
$ seqkit seq --name --only-id read_1.fq.gz read_2.fq.gz \
| sort | uniq -d > id.txt
$ # number of IDs
wc -l id.txt
8090 id.txt
## 然后用id.txt来提取reads
$ seqkit grep --pattern-file id.txt read_1.fq.gz -o read_1.f.fq.gz
$ seqkit grep --pattern-file id.txt read_2.fq.gz -o read_2.f.fq.gz
## 用md5sum来检查两个文件的IDs是否是相同的
$ seqkit seq --name --only-id read_1.f.fq.gz > read_1.f.fq.gz.id.txt
$ seqkit seq --name --only-id read_2.f.fq.gz > read_2.f.fq.gz.id.txt
$ md5sum read_*.f.fq.gz.id.txt
537c57cfdc3923bb94a3dc31a0c3b02a read_1.f.fq.gz.id.txt
537c57cfdc3923bb94a3dc31a0c3b02a read_2.f.fq.gz.id.txt
##sort 一下
$ gzip -d -c read_1.f.fq.gz \
| seqkit fx2tab \
| sort -k1,1 -T . \
| seqkit tab2fx \
| gzip -c > read_1.f.sorted.fq.gz
$ gzip -d -c read_2.f.fq.gz \
| seqkit fx2tab \
| sort -k1,1 -T . \
| seqkit tab2fx \
| gzip -c > read_2.f.sorted.fq.gz
将两个fasta文件连接成一个
$ cat 1.fa
>seq1
aaaaa
>seq2
ccccc
>seq3
ggggg
$ cat 2.fa
>seq3
TTTTT
>seq2
GGGGG
>seq1
CCCCC
一行命令
$ seqkit concat 1.fa 2.fa
>seq1
aaaaaCCCCC
>seq2
cccccGGGGG
>seq3
gggggTTTTT