生信单行脚本

本文总结了,生物信息处理过程中常见一些工具和单行命令等。

awk和sed基础

提取文件中的2, 4, and 5 列:

输出第五列等于abc123的行:

输出第五列不是abc123的行:

输出第七列以字母a-f开头的行:

输出第七列不是以字母a-f开头的行:

计算第二列不重复的值保存在哈希arr中 (一个值只保存一次):

输出第三列的值比第五列大的行:

计算文件中第一列的累加值,输出最后的结果:

计算第二列的平均值:

用bar替换文件中所有的foo:

消除行开头空和格制表符:

消除行结尾的空格和制表符:

消除行中开头和结尾的空格和制表符:

删除空行:

删除包含‘EndOfUsefulData’的行及其后所有的行:

生信sed,awk单行应用

Returns all lines on Chr 1 between 1MB and 2MB in file.txt. (assumes) chromosome in column 1 and position in column 3 (this same concept can be used to return only variants that above specific allele frequencies):

输出Chr为1在1M和2M之间的所有行。(假设)染色体在第一列,位点在第三列(基于同样的假设可以用来返回类似特定等位基因频率的变异)

Basic sequence statistics. Print total number of reads, total number unique reads, percentage of unique reads, most abundant sequence, its frequency, and percentage of total in file.fq: 基本序列统计。输出总的reads数,不重复的reads总数,不重复reads百分比,最大冗余的序列及其频度以及总占比百分数。

cat myfile.fq | awk ‘((NR-2)%4==0){read=$1;total++;count[read]++}END{for(read in count){if(!max||count[read]>max) {max=count[read];maxRead=read};if(count[read]==1){unique++}};print total,unique,unique*100/total,maxRead,count[maxRead],count[maxRead]*100/total}’

转换.bam为.fastq:

Keep only top bit scores in blast hits (best bit score only): 只取blast采样中的顶级位点的分数(最高的位点分)

Keep only top bit scores in blast hits (5 less than the top): 只取blast采样中的顶级位点的分数(比顶级少于5的)

分割多序列FASTA文件为单序列FASTA文件

输出fasta文件中的每条序列的序列名称和长度

转化FASTQ文件为FASTA:

从第二行开始每四行取值(从FASTQ文件提取序列)。

输出中剔除第一行:

输出20-80行:

计算二,三行列的和并追加到每行后输出

计算fastq文件平均reads的长度

转化VSF文件为BED文件

sort, uniq, cut等杂项

输出带行号的内容:

去重复行计数

找到两文件都有的行(假设两个文件都是无重复行,重定向执行‘wd -l’计算同样行的行数)

对文件按照第九列数字顺序排序(g按照常规数值,k列)

找到第二列出现最多的字符串

从文件中随机取10行

输出所有三个所可能的DNA序列

Untangle an interleaved paired-end FASTQ file. If a FASTQ file has paired-end reads intermingled, and you want to separate them into separate /1 and /2 files, and assuming the /1 reads precede the /2 reads:

解开一列交错paired-end fastq文件。如果fastq文件有乱序paired-end reads,你想将其分离成单独的/1,/2的文件保存,这里假设/1 reads 在/2 前面:

Take a fasta file with a bunch of short scaffolds, e.g., labeled >Scaffold12345, remove them, and write a new fasta without them:

将一个fasta文件转成一系列短的scaffolds。比如,标签 “>Scaffold12345″,然后移出他们,保存一个去掉他们的新文件:

Display hidden control characters:

显示一个隐藏的控制字符:

find, xargs,和GNU parallel

通过 https://www.gnu.org/software/parallel/. 载 GNU parallel

搜索文件夹及其子目录中名称为 .bam 文件(目录也算):

删除上面搜到的文件列表(不可逆的危险操作,谨慎使用!删除之前请自习确认)

将所有.txt 文件修改为.bak(例如在对*.txt做操作之前用于文件备份)

Chastity filter raw Illumina data (grep reads containing :N:, append (-A) the three lines after the match containing the sequence and quality info, and write a new filtered fastq file):

对Illumina数据做Chastity过滤(grep 查询 包含:N:,用(-A)选项第三列信息附加在匹配的包含一个序列质量信息后,并保存为一个新的fasta文件)

通过parallel并行运行12个FASTQC任务

通过parallel给bam做索引,通过–dry-run打印测试这些命令,实际上并未做执行。

seqtk

Seqtk项目托管地址https://github.com/lh3/seqtk。Seqtk是一个快捷轻量的处理FASTA和FASTQ格式基因序列的工具。他可以是先FASTA和FASTQ无缝处理和转化,同时支持gzip格式的压缩文件。

把FASTQ转化为FASTA:

转化ILLUMINA 1.3+ 格式FASTQ为FASTA,并且以小于20的mask bases获得小写字母(第一命令行)或者到N(第二)。 seqtk seq -aQ64 -q20 in.fq > out.fa seqtk seq -aQ64 -q20 -n N in.fq > out.fa

折叠长FASTA/Q行,并且去除其注释:

转化多行FASTQ到四行FASTQ:

反转FASTA/Q序列:

用序列文件中的名称(比如name.1st)提取序列,一个虚列名一行:

利用序列文件中的”reg.bed“r信息提取地理信息的序列:

编码‘reg.bed’信息为小写

从两个大的paired FASTQ文件提取10000个read pairs(记得用同样的随机种子保持 paire)

利用Phred公式从两头修剪低质量bases:

从左端修剪5bp,从右端修剪10bp的。

GFF3 Annotations

输出GFF3文件中标注的所有的序列

检测GFF3文件中标注的所有性状类型。

检测GFF3文件中标注的基因数量。

从GFF3文件中提取所有的基因ID

输出GFF3文件每个基因的长度

FASTA头列转化为GFF格式(假设头的长度,附加在”_length“ ,和Velvet assembled transcripts))

有用的别名(.bashrc)

提示符修改为[email protected]:/full/path/cwd/:$ 形式

export PS1=”\[email protected]\h:\w\\$ ”

避免反复敲诸如cd ../../..的命令(也可以用[autojump](https://github.com/joelthelion/autojump),让你在飞速的转换目录

向前和向后浏览

覆盖文件时候,先确认

我最喜欢的”ls“别名

对cut空格和逗号,分割文件

解压缩tar包

或者可以用更普遍的‘extract’函数

# 源于ABSG(Advanced Bash Scripting Guide)中 Mendel Cooper的建议

使用别名”mcd”创建一个目录,并且cd到该目录

跳转到上级目录,并且列出其内容

一个好看的grep

刷新你的.bashrc

编辑你的.bashrc

常用错误别称

使用 pandoc转化markdown文档为PDF格式:

# 用法: mdpdf document.md document.md.pdf

对当前目录搜索关键词(ft “mytext” *.txt):

Etc

重复运行上一条命令:

sudo !!

敲出了部分命令,删除这些输入,查你忘记的明亮,拉回命令,继续输入(<CTRL+u>删除光标之前的输入,<CTRL+y>恢复上个C-U删除字符)

跳到一个目录,执行命令,然后返回当前目录(()的用法)

记时秒表 (输入Enter or ctrl-d 停止):

把上次执行的命令生成一个脚本

重用上次命令的所有参数

列出或者删除一个目录中所有不匹配的特定后缀的文件(例如,列出所有不是压缩的文件,删除所有不以.foo和.bar后缀的文件)

利用上次的命令,但是不需要他的的参数(重新输入参数):

激活一个快捷的编辑器,输入,编辑长的,复杂,巧妙的命令:

输出一个特定的行(比如 42行)

终结一个冻结的ssh session(会车换行,敲~键,在敲下.键)

利用grep去除文件的空行,结果保存到新文件

查找大文件(例如,大于500M的)

利用截取列(例如,一个tab分割文件的第五个域)

查找包含特定字符的文件(-l 只输出文件名, -i 忽略大小写 -r 遍历子目录)

发表评论

邮箱地址不会被公开。 必填项已用*标注