[基因组工具]seqkit的使用

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
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,732评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,496评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,264评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,807评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,806评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,675评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,029评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,683评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 41,704评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,666评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,773评论 1 332
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,413评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,016评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,978评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,204评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,083评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,503评论 2 343

推荐阅读更多精彩内容