【linux】单行命令-第2部分:Bioinformatics one-liners

单行命令,化繁为简,重剑无锋,大巧不工。接下来学习生信实用单行命令。

基础awk和sed命令

sed: stream editor流编辑器

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

awk '{print $2,$4,$5}' input.txt

输出文件第5列中等于abc123的行:

awk ` $5 == "abc123" ' input.txt

输出第5列不等于abc123的行:

awk '$5 != "abc123" ' input.txt

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

awk '$7 ~ /^[a-f]/' input.txt

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

awk '$7 !~ /^[a-f]/' input.txt

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

awk '!arr[$2]++' input.txt

输出第3列值比第5列大的行:

awk '$3>$5' input.txt

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

awk '{sum+=$1} END {print sum}'  input.txt

计算第2列的平均值:

awk '{x+=$2} END {print x/NR}' input.txt

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

sed 's/foo/bar/g' input.txt

消除行首空格或制表符:

sed 's^[ \t]*//' input.txt

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

sed 's/[ \t]*$//' input.txt

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

sed 's/^[ \t]*//;s/[ \t]*$//' input.txt

删除空行:

sed '/^$/d' input.txt

删除包含“EndOfUsefulData”的行及其后所有的行:

sed -n '/EndOfUsefulData/,$!p' input.txt

生信单行sed&awk

提取file.txt 文件Chr1中1MB和2MB的片段之间的行信息,假设第1列是Chr信息,第3列是位置信息:

# bed 文件
cat file.bed | awk '$1 == "1" ' | awk '$3>=999999' | awk '$3<=1999999'
# gff文件
cat file.gff | awk '$1 == "1" '| awk '$5 >=1000000' | awk '$5<=2000000'

统计fastq文件的一些基本信息,包括reads数、唯一的reads数、唯一reads数的比例、出现频率最多的序列及频数和所占比例:(拿到rawdata来统计时很实用)

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文件:

samtools view file.bam | awk 'BEGIN {FS='\t'} {print '@' $1 '\n' $10 '\n+\n' $11}' > file.fq

输出blast结果中最高得分的结果:

awk '{if (!x[$1]++) {print $0; bitscore=($14-1)} else {if($14>bitscore) print $0} }' blastout.txt

将含多条序列fasta文件分割为多个fasta(每条一个文件):

awk '/^>/{s=++d".fa"} {print > s}' multi.fa

输出fasta文件中每条序列的序列名及其长度:

cat file.fa | awk `$0 ~ ">"{print c;c=0; printf substr($0,2,100) "\t"; } $0 !~ ">" {c+=length($0);} END { print c; }'

将fastq文件转为fasta文件:

sed -n '1~4s/^@/>/p;2~4p' file.fq > file.fa

从第2行开始,每四行取值(从FASTQ文件提取序列):

sed -n '2~4p' file.fq

输出中剔出第1行:

awk 'NR>1' input.txt

输出第20-80行:

awk 'NR>=20&&NR<=80' input.txt

计算第2、3行列的和,并追加到每行后输出:

awk '{print $0, $2+$3}' input.txt

计算fastq文件评价read的长度

awk  'NR%4==2(sumi+=length($0)}END{print sum(NR/4)}' input.fastq

转化VSF文件为BED文件:

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

sort, uniq,cut

输出带行号的内容:

cat -n file.txt

去除重复行并计数:

cat file.txt | sort -u | wc -l

找到两文件都有的行(如果两文件都是无重复行,重定向执行“wd -l”计算同样行的行数)

sort file1 file2 | uniq -d

# 安全的方法
sort -u file1 > a
sort -u file2 > b
sort a b | uniq -d 

# 用comm的方法
comm -12 file1 file2

文件按照第9列数字排序(g按照常规数值,k列):

sort -gk9 file.txt

找到第2列出现最多的字符:

cut -f2 file.txt | sort |uniq -c | sort -k1nr | head

从文件中随机选取10行:

shuf file.txt | head -n 10

输出所有的3mer DNA组合:

echo {A,C,T,G}{A,C,T,G}{A,C,T,G}

将合并的paired-end fastq文件拆分为-1和-2两个文件:(这里加上/1在/2前面)

cat interleaved.fq | paste - - - - - - - - | tee > (cut -f 1-4 | tr "\t" "\n" > deinterleaved_1.fq) |cut -f 5-8 | tr "\t" "\n" >deinterleaved_2.fq

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

samtools faidx genome.fa && grep -v Scaffold genome.fa.fai | cut -f1 | xargs -n1 samtools faidx genome.fa > geome.noscaffolds.fa

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

python -c "f = open ('file.txt' , 'r' ); f.seek(0); file = f.readlines(); print file"

find, xargs, GNU parallel

xargs:“eXtended ARGuments”的缩写,主要与find,echo和cp等命令结合使用。
find:在指定目录下查找。

GNU parallel 下载地址:https://www.gnu.org/software/parallel/
搜素文件夹及其子目录中名称为.bam的文件(目录也包括在内):

find . -name "*.bam"

删除所有的bam文件:(谨慎操作!!!)

find . -name "*.bam" | xargs rm

将所有.txt文件重命名为.bak:(如在对*.txt做操作前用于文件备份)

find . -name "*.txt" | sed "s/\.txt$//" | xargs -i echo mv {}.txt {}.bak | sh

同时运行12个FASTQC文件:

find *fq | parallel -j 12 "fastqc {} --outdir . "

将bam文件建索引(进输出命令,并不进运行程序):

find *.bam | parallel --dry-run 'samtools index {}'

seqktk

Seqtk专为FASTA和FASTQ而生,能快速处理FASTA或FASTQ格式序列,也可读取gzip压缩过的FASTA和FASTQ文件。下载地址:https://github.com/lh3/seqtk

将fastq转为fasta格式:

seqtk seq -a in.fq.gz > out.fa

将fastq文件中的质量值低于20的序列屏蔽掉并转为fasta格式:

# 将序列屏蔽为小写
seqtk seq -aQ64 -q20 in.fq > out.fa
# 将序列屏蔽为N
seqtk seq -aQ64 -q20 -n N in.fq > out.fa

将fasta和fastq文件格式化为每行60个字符的多行序列并去除注释信息:

seqtk seq -Cl60 in.fa > out.fa

将多行的fastq文件转为4行的fastq文件:

seqtk seq -l0 in.fq > out.fq

生成fastq或fasta的反向互补序列:

seqtk seq -r in.fq > out.fq

根据name.lst(每行一个序列名)中的序列名提取序列:

seqtk subseq in.fq name.lst > out.fq

提取reg.bed文件中所含区域的序列:

seqtk subseq in.fa reg.bed > out.fa

将reg.bed文件中所含区域的序列屏蔽为小写:

seqtk seq -M reg.bed in.fa > out.fa

使用Phred算法从两端修剪低质量的碱基:

seqtk trimfq in.fq > out.fq

从每条read的左端修剪5bp,从右端修剪10bp:

seqtk trimfq -b 5 -e 10 in.fa > out.fa

将合并的paired-end fastq文件拆分为-1和-2 两个fastq文件:

seqtk seq -l0 -1 interleaved.fq > deinterleaved_1.fq
seqtk seq -l0 -2 interleaved.fq > deinterleaved_2.fq

GFF3注释文件

输出GFF3文件所有注释的序列:

cut -s -f 1,9 yourannots.gff3 | grep $ '\t' | cut -f 1 | sort | uniq

统计GFF3文件的基因数量:

grep -c $'\tgene\t'  yourannots.gff3

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

grep $'\tgene\t' yourannots.gff3 | perl -ne 'ID=([^;]+)/ and  printf("%s\n", $1)'

统计GFF3文件每个基因的长度

grep $'\tgene\t' yourannots.gff3 | cut -s -f 4,5 | perl -ne '@v = split(/\t); printf("%d\n", $v[1]-$v[0]+1)'

参考链接:

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

推荐阅读更多精彩内容