【生信工具】seqkit-fa&fq文件处理利器

作为生信入门训练,我们常用perl python等脚本语言实现对基因组文件的处理,练习常规的文本文件处理。最近再做单细胞方面的技术开发,需要自己手动合成一个嵌合两种物种的基因组,且要配套对应的gtf文件。用于评估单细胞技术的双胞率。
于是把seqkit这个由重庆第三军医大学Wei Shen 使用 go 语言开发用于处理 fa 和 fq 文件的利器拿出来温习一下,发现我的需求基本都可以满足。同等于处理bed文件的万能工具bedtools,感慨一套处理工具软件包可以取代一批生信工程师,这话确实如此。

主要功能
一、序列操作:

1.取反向序列

seqkit seq test.fa -r > test_re.fa

2.取互补序列

seq test.fa -p > test_com.fa

3.取反向互补序列

seqkit seq test.fa -r -p > test_re_com.fa

4.DNA序列转换为RNA序列

seqkit seq test.fa --nda2rna > test_rna.fa

5.RNA序列转换为DNA序列

seqkit seq test.fa rna2dna > test_dna.fa

6.将序列以小写字母的形式输出

seqkit seq test.fa -l > test_lower.fa

7.将序列以大写字母的形式输出

seqkit seq test.fa -u > test_upper.fa

8.指定每行序列的输出长度(为0的话,代表为一整行,默认的输出 长度是60个碱基)

seqkit seq test.fa -w 50 > test_50.fa (指定序列的长度为50)

9.将多行序列转换为一行序列

seqkit seq test.fa -w 0 > test_w.fa

10.只输出序列

seqkit seq test.fa -s -w 0 > test_seq.fa

11.将只输出的序列的,指定每行输出的碱基数

seqkit seq test_seq.fa -s -w 40 > test_seq40.fa

12输出参数-o
seqkit seq test.fa -s -w 20 -o test_20.fa

二、Fasta/q之间以及与tab格式互换

1.将fataq文件转化为fasta格式.

seqkit seq fq2fa test.fq -o test.fa

2.将fasta格式转化为tab格式

seqkit fx2tab test.fa > test_tab.fa (没有seq参数)

三、序列信息统计

1.序列碱基含量

seqkit fx2tab -l -g -n -i -H test.fa (这些参数组合起来比较好看)

2.序列长度的整体分布统计

seqkit stat test.fa
seqkit grep [flags]    #这部分内容用的比较广泛

参数:

-n, --by-name

匹配整个序列的名字,包含deion部分,而不是序列id。

-s, --by-seq

匹配序列

-d, --degenerate

pattern/motif 包含简并碱基

-i, --ignore-case

忽略大小写

-v, --invert-match

输出不匹配此模式的内容

-p,

匹配模式,支持连续写多个模式,匹配任一模式即输出。如-p ^ATG -p TAA$。注意该功能仅能正向匹配,不能实现对互补链匹配。

-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 list 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##匹配限定到某区域

五、motif定位

对grep的拓展,可以正反链同时匹配,输出匹配的位置。

seqkit locate [flags]

参数:

-d, --degenerate

pattern/motif contains degenerate base

-i, --ignore-case

ignore case

-P, --only-positive-strand

only search at positive strand

-p, --pattern value

search pattern/motif

-f, --pattern-file string

pattern/motif file (FASTA format)

举例:

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

六、多个序列文件比较寻找相同的序列或者ID相同的序列

seqkit common [flags]

参数:

-n, --by-name

匹配整个序列的名字,包含deion部分,而不是序列id

-s, --by-seq

match by sequence

-i, --ignore-case

ignore case

-m, --md5

use MD5 reduce memory usage

举例:

1、By ID (default,>后面,空格之前的名字)输出ID名字相同的。

seqkit common test1.fa test2.fa -o common.fasta

2、By full name(整个序列的名字,包含deion部分)。输出序列名字相同的。

seqkit common test1.fa test2.fa -n -o common.fasta

3、输出要比较的文件中序列相同的序列

seqkit common test1.fa test2.fa -s -i -o common.fasta

4、输出要比较的文件中序列相同的序列 (for large sequences)

seqkit common test1.fa test2.fa -s -i -o common.fasta --md5

七、提取部分序列

如随机抽取10000条FASTQ序列做NT污染评估。同时他也可以对FASTA序列提取

seqkit sample [flags]

参数:

-n, --number int

sample by number (result may not exactly match)

-p, --proportion float

sample by proportion(按比例提)

-s, --rand-seed int

rand seed for shuffle (default 11)

-2, --two-pass

2-pass modelower memory

举例:随机抽取序列

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

八、排序输出命令

seqkit sort [flags]

参数:

-l, --by-length

按照序列长度排序

-n, --by-name

by full name

-s, --by-seq

按照序列排序

-i, --ignore-case

按序列排序时忽略大小写

-r, --reverse

反向排序

-2, --two-pass

对于FASTA序列排序可以减少内存

举例:

seqkit sort -ltest.fa

九、文件切割

seqkit split [flags]

参数:

-i, --by-id

split squences according to sequence ID

-p, --by-part int

将一个文件分割成N 份

-s, --by-size int

将一个文件按照N 条序列一个文件进行分割

-O, --out-dir string

output directory (default value is infile.split)

-2, --two-pass

two-pass mode to lower memory usage(only FAST)

举例:

seqkit split hairpin.fa.gz -p 4
参考材料:

1.https://github.com/shenwei356/seqkit
2.https://www.cnblogs.com/huangyinger/p/10421805.html

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

推荐阅读更多精彩内容