分类目录归档:shell

生信单行脚本

单行脚本处理生信息常见脚本

Basic awk & sed

[返回顶部]

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

输出第五列等于abc123的行:

输出第五列不是abc123的行:

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

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

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

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

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

计算第二列的平均值:

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

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

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

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

删除空行:

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

awk & sed for bioinformatics

生信单行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百分比,最大冗余的序列及其频度以及总占比百分数。

转换.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文件

sed -e ‘s/chr//’ file.vcf | awk ‘{OFS=”\t”; if (!/^#/){print $1,$2-1,$2,$4″/”$5,”+”}}’

sort, uniq, cut, etc.

[返回开头]

输出带行号的内容:

去重复行计数

找到两文件都有的行(假设两个文件都是无重复行,重定向执行‘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, and 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打印测试这些命令,实际上并未做执行。 find *.bam | parallel –dry-run ‘samtools index {}’

seqtk

[back to top]

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

把FASTQ转化为FASTA:

Convert ILLUMINA 1.3+ FASTQ to FASTA and mask bases with quality lower than 20 to lowercases (the 1st command line) or to N (the 2nd):

转化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

Fold long FASTA/Q lines and remove FASTA/Q comments:

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

Convert multi-line FASTQ to 4-line FASTQ: 转化多行FASTQ到四行FASTQ:

Reverse complement FASTA/Q: 反转FASTA/Q序列:

Extract sequences with names in file name.lst, one sequence name per line: 用序列文件中的名称(比如name.1st)提取序列,一个虚列名一行:

Extract sequences in regions contained in file reg.bed: 利用序列文件中的”reg.bed“r信息提取地理信息的序列:

Mask regions in reg.bed to lowercases: 编码‘reg.bed’地理信息到小写

Subsample 10000 read pairs from two large paired FASTQ files (remember to use the same random seed to keep pairing): 从两个大的paired FASTQ文件提取10000个read pairs(记得用同样的随机种子保持 paire)

Trim low-quality bases from both ends using the Phred algorithm: 利用Phred公式从两头修剪低质量bases:

Trim 5bp from the left end of each read and 10bp from the right end: 从左端修剪5bp,从右端修剪10bp的。

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 read在/2 read的前面): seqtk seq -l0 -1 interleaved.fq > deinterleaved_1.fq seqtk seq -l0 -2 interleaved.fq > deinterleaved_2.fq

GFF3 Annotations

[back to top]

Print all sequences annotated in a GFF3 file. 输出GFF3文件中标注的所有的序列

Determine all feature types annotated in a GFF3 file. 检测GFF3文件中标注的所有性状类型。

Determine the number of genes annotated in a GFF3 file. 检测GFF3文件中标注的基因数量。

Extract all gene IDs from a GFF3 file. 从GFF3文件中提取所有的基因ID

Print length of each gene in a GFF3 file. 输出GFF3文件每个基因的长度

FASTA header lines to GFF format (assuming the length is in the header as an appended “_length” as in Velvet assembled transcripts): FASTA头列转化为GFF格式(假设头的长度,附加在”_length“ ,和Velvet assembled transcripts))

Other generally useful aliases for your .bashrc

有用的别名(.bashrc)

[back to top]

提示符修改为user@hostname:/full/path/cwd/:$形式

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

向前和向后浏览

覆盖文件时候,先确认

我最喜欢的”ls“别名

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

解压缩tar包

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

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

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

一个好看的grep

刷新你的.bashrc

编辑你的.bashrc

常用错误别称

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

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

Etc

[返回]

重复运行上一条命令:

列出最近最常用的命令行参数(通常是文件)

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

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

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

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

重用上次命令的所有参数

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

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

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

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

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

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

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

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

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

nginx日志自动切割以清理

按日自动切割nginx日志,并删除30(可自定义)天以上的日志。
脚本保存为 /nginx/logs/ngx_logcut.sh 并加入crontab中

iptables和ipset批量限制非法ip源

一、获取非法ip源:

扫描登录失败的日志:
less secure*|grep “Failed password for root from”  >/tmp/gj.list
对扫描日志进行分析,非法尝试失败超过50的ip列出来,保存为gongji.ip
cat /tmp/gj.list|perl -anle  ‘print $F[10]’|sort |uniq -c|sort -n |perl -anle ‘print $F[1] if $F[0]>50 ‘>/tmp/gongji.ip

继续阅读