Super useful!! #实用的单行命令合集

周末闲着无聊在推特上闲逛，刷着刷着发现了Ming Tang大神的一条推文。他将他日常经常使用的单行命令整理到其Github中，然后分享给各位推友。在点开它的网站后，我犹如发现了一堆宝藏。里面很多单行命令行非常实用，快速简单地解决了我日常所遇到的生物信息学数据处理问题(没有这些单行命令，我往往要通过自己写脚本或者借用其它工具来解决日常的生物信息学数据处理)。在取得Ming Tang大神的同意后，该推文将和大家整理分享这一系列实用的单行命令。

原文链接：
GitHub 地址： https://github.com/<wbr>crazyhottommy/bioinformatics-<wbr>one-liners

Fastq/Fasta文件的处理

统计fastq文件中的read的长度与不同长度的分布：

zcat file.fastq.gz | awk 'NR%4 == 2 {lengths[length($0)]++} END {for (l in lengths) {print l, lengths[l]}}'

fastq转换为fasta格式:

zcat file.fastq.gz | paste - - - - | perl -ane 'print ">$F[0]\n$F[2]\n";' | gzip -c > file.fasta.gz

统计fastq文件中reads的数目:

cat file.fq | echo $((`wc -l`/4))

随机抽取fastq文件的一部分，下面使用提取0.01%的reads作为例子(在不知道该方法之前，我一般使用Heng Li的seqtk sample来达到该目的):

cat file.fq | paste - - - - | awk 'BEGIN{srand(1234)}{if(rand() < 0.01) print $0}' | tr '\t' '\n' > out.fq

从fastq文件中提取含有某个碱基排列序列的read，并且让其保持fastq的格式(含有四行的信息:

### 方法一
grep -A 2 -B 1 'AAGTTGATAACGGACTAGCCTTATTTT' file.fq | sed '/^--$/d' > out.fq

### 方法二
zcat reads.fq.gz \
| paste - - - - \
| awk -v FS="\t" -v OFS="\n" '$2 ~ "AAGTTGATAACGGACTAGCCTTATTTT" {print $1, $2, $3, $4}' \
| gzip > filtered.fq.gz

分选fastq文件，有些fastq文件是同时包含了R1,R2的序列信息，需要将它们分离出来:

cat file.fq | paste - - - - - - - - | tee >(cut -f1-4 | tr '\t'  
'\n' > out1.fq) | cut -f5-8 | tr '\t' '\n' > out2.fq

将一个含有多个染色体或者contigs的fasta文件，按照染色体或者contigs分离成多个只含单一序列的fasta文件:

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

将fasta文件中多行的序列变成单行:

### 方法一
cat file.fasta | awk '/^>/{if(N>0) printf("\n"); ++N; printf("%s\t",$0);next;} {printf("%s",$0);}END{printf("\n");}'

### 方法二
awk 'BEGIN{RS=">"}NR>1{sub("\n","\t"); gsub("\n",""); print RS$0}' file.fa

反过来将单行的fasta文件，变成多行，下面例子对每60个碱基分一行:

awk -v FS= '/^>/{print;next}{for (i=0;i<=NF/60;i++) {for (j=1;j<=60;j++) printf "%s", $(i*60 +j); print ""}}' file

fold -w 60 file

统计mulitfasta中，每个染色体/contigs的长度:

awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen = seqlen +length($0)}END{print seqlen}' file.fa

对fasta文件中的header中添加对应的信息:

### 目前你有的fasta文件
cat myfasta.txt 
>Blap_contig79
MSTDVDAKTRSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>Bluc_contig23663
MSTNVDAKARSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>Blap_contig7988
MSTDVDAKTRSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>Bluc_contig1223663
MSTNVDAKARSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI

###需要添加的信息
cat my_info.txt 
info1
info2
info3
info4

### 添加信息到对应的fasta header中
paste <(cat my_info.txt) <(cat myfasta.txt| paste - - | cut -c2-) | awk '{printf(">%s_%s\n%s\n",$1,$2,$3);}'
>info1_Blap_contig79
MSTDVDAKTRSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>info2_Bluc_contig23663
MSTNVDAKARSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>info3_Blap_contig7988
MSTDVDAKTRSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>info4_Bluc_contig1223663
MSTNVDAKARSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI

简化fasta文件的表头,比如将>7 dna:chromosome chromosome:GRCh37:7:1:159138663:1转化为>7

cat Homo_sapiens_assembly19.fasta | gawk '/^>/ { b=gensub(" dna:.+", "", "g", $0); print b; next} {print}' > Homo_sapiens_assembly19_reheader.fasta

BAM/VCF/BED文件的处理

将BAM文件转化成bed格式

samtools view file.bam | perl -F'\t' -ane '$strand=($F[1]&16)?"-":"+";$length=1;$tmp=$F[5];$tmp =~ s/(\d+)[MD]/$length+=$1/eg;print "$F[2]\t$F[3]\t".($F[3]+$length)."\t$F[0]\t0\t$strand\n";' > file.bed

将BAM文件转化成wig格式

samtools mpileup -BQ0 file.sorted.bam | perl -pe '($c, $start, undef, $depth) = split;if ($c ne $lastC || $start != $lastStart+1) {print "fixedStep chrom=$c start=$start step=1 span=1\n";}$_ = $depth."\n";($lastC, $lastStart) = ($c, $start);' | gzip -c > file.wig.gz

由于samtools mpileup只能使用1个CPU，按照染色体，并使用xargs去并行call variant能够极大地提高其变异发现的效率:

samtools view -H yourFile.bam | grep "\@SQ" | sed 's/^.*SN://g' | cut -f 1 | xargs -I {} -n 1 -P 24 sh -c "samtools mpileup -BQ0 -d 100000 -uf yourGenome.fa -r {} yourFile.bam | bcftools view -vcg - > tmp.{}.vcf"

另外在分染色体call完变异后，也可以将不同染色体的变异合并起来:

samtools view -H yourFile.bam | grep "\@SQ" | sed 's/^.*SN://g' | cut -f 1 | perl -ane 'system("cat tmp.$F[0].bcf >> yourFile.vcf");'

寻找当前目录下的bam文件，并使用5个CPUs将其复制到其它新的目录中(这个命令在转移文件时，非常实用):

find . -name "*bam" | xargs -P5 -I{} rsync -av {} dest_dir

根据vcf的表头文件来sort VCF文件

cat my.vcf | awk '$0~"^#" { print $0; next } { print $0 | "sort -k1,1V -k2,2n" }'

将含有PASS的变异从VCF中提取出来(经过一些工具过滤后，例如GATK，如果符合过滤的阀值，在VCF文件中会有一行PASS的字符):

cat my.vcf | awk -F '\t' '{if($0 ~ /\#/) print; else if($7 == "PASS") print}' > my_PASS.vcf

将bed文件按照不同的染色体分开(同样的道理可以应用到gff文件中):

### 方法一
cat nexterarapidcapture_exome_targetedregions_v1.2.bed | sort -k1,1 -k2,2n | sed 's/^chr//' | awk '{close(f);f=$1}{print > f".bed"}'

### 方法二
awk '{print $0 >> $1".bed"}' example.bed

处理其它文件/序列

将序列转化其互补的序列，在设计primers很有用:

echo 'ATTGCTATGCTNNNT' | rev | tr 'ACTG' 'TGAC'

重命名一系列文件:

for file in *gz
do zcat $file > ${file/bed.gz/bed}

将看不到的字符显示出来

### 方法一
cat my_file | sed -n 'l'
### 方法二
cat -A

复制大量的文件，使用rsync比copy速度更快

## 将所有文件从一个目录复制到另一个目录
rsync -av from_dir/ to_dir
##将所有文件从一个目录复制到另一个目录，但是已经复制好的文件将会被跳过:
rsync -avhP /from/dir /to/dir

显示当前目录下所有的子目录所占空间的大小

du -h --max-depth=1

查看当前内存的使用情况

free -mg

根据第一和第二列打印出特殊唯一的行:

awk '!a[$1,$2]++' input_file
### 或者
sort -u -k1,2 file

删除空白行:

sed /^$/d

删除最后一行

sed $d

使用!$,使用最后一个命令的字符

echo hello, world; echo !$

将最后执行的命令保存为脚本:

echo "!!" > foo.sh

重复使用刚刚执行命令的参数:

!*

统计一个tsv文件中有多少列:

cat file.tsv | head -1 | tr "\t" "\n" | wc -l  

##csvkit工具中的csvcut
csvcut -n -t file.

## 另一些可行的方法
less files.tsv | head -1| tr "\t" "\n" | nl
awk -F "\t" 'NR == 1 {print NF}' file.tsv
awk '{print NF; exit}'

创建一个新的文件夹，并且马上cd进去

mkdir blah && cd $_

切割所需要的行，根据另一个文件提供的行名称(Thanks God,我一直想找现成的办法去处理这个问题,非常合适处理大型的文件，比R占更少的内存):

#! /bin/bash

set -e
set -u
set -o pipefail

#### Author: Ming Tang (Tommy)
#### Date 09/29/2016
#### I got the idea from this stackOverflow post http://stackoverflow.com/questions/11098189/awk-extract-columns-from-file-based-on-header-selected-from-2nd-file

# show help
show_help(){
cat << EOF
  This is a wrapper extracting columns of a (big) dataframe based on a list of column names in another
  file. The column names must be one per line. The output will be stdout. For small files < 2G, one 
  can load it into R and do it easily, but when the file is big > 10G. R is quite cubersome. 
  Using unix commands on the other hand is better because files do not have to be loaded into memory at once.
  e.g. subset a 26G size file for 700 columns takes around 30 mins. Memory footage is very low ~4MB.

  usage: ${0##*/} -f < a dataframe  > -c < colNames> -d <delimiter of the file>
        -h display this help and exit.
        -f the file you want to extract columns from. must contain a header with column names.
        -c a file with the one column name per line.
        -d delimiter of the dataframe: , or \t. default is tab.  
        
        e.g. 
        
        for tsv file:
            ${0##*/} -f mydata.tsv -c colnames.txt -d $'\t' or simply ommit the -d, default is tab.
        
        for csv file: Note you have to specify -d , if your file is csv, otherwise all columns will be cut out.
            ${0##*/} -f mydata.csv -c colnames.txt -d ,
        
EOF
}

## if there are no arguments provided, show help
if [[ $# == 0 ]]; then show_help; exit 1; fi

while getopts ":hf:c:d:" opt; do
  case "$opt" in
    h) show_help;exit 0;;
    f) File2extract=$OPTARG;;
    c) colNames=$OPTARG;;
    d) delim=$OPTARG;;
    '?') echo "Invalid option $OPTARG"; show_help >&2; exit 1;;
  esac
done
    

## set up the default delimiter to be tab, Note the way I specify tab 

delim=${delim:-$'\t'}

## get the number of columns in the data frame that match the column names in the colNames file.
## change the output to 2,5,6,22,... and get rid of the last comma  so cut -f can be used
 
cols=$(head -1 "${File2extract}" | tr "${delim}" "\n" | grep -nf "${colNames}" | sed 's/:.*$//' | tr "\n" "," | sed 's/,$//')

## cut out the columns 
cut -d"${delim}" -f"${cols}" "${File2extract}"

相同的问题，使用csvtk就简单很多了:

csvtk cut -t -f $(paste -s -d , list.txt) data.tsv

在每行开始的地方增加或者删除chr，有些软件不支持chr字符，需要使用数字来表示染色体:

# 增加 chr
sed 's/^/chr/' my.bed

# 或者
awk 'BEGIN {OFS = "\t"} {$1="chr"$1; print}'

# 移除 chr
sed 's/^chr//' my.bed

检查 tsv文件中行与列的数目的是否一致:

awk '{print NF}' test.tsv | sort -nu | head -n 1

只替代其中一列中的某个pattern

## 替代第五列
awk '{gsub(pattern,replace,$5)}1' in.file

## 使用强大的csvtk
csvtk replace -f 5 -p pattern -r replacement

反转文件中的其中一列

###反转第三列然后将其结果放到第五列中
awk -v OFS="\t" '{"echo "$3 "| rev" | getline $5}{print $0}' 
###只反转第二列
perl -lane 'BEGIN{$,="\t"}{$rev=reverse $F[2];print $F[0],$F[1],$rev,$F[3]}

从gtf文件中提取promoter的位置，例子提取基因上下游5kbp的序列

zcat gencode.v29lift37.annotation.gtf.gz | awk '$3=="gene" {print $0}' | grep protein_coding | awk -v OFS="\t" '{if ($7=="+") {print $1, $4, $4+5000} else {print $1, $5-5000, $5}}' > promoters.bed

但是值得注意的是，某些基因位于染色体的末端，提取上下游5kbp可能超出了这一点，这种情况下最好请使用bedtools slop根据基因组大小文件来避免这种情况。

## 使用fetchChromSizes提取染色体大小，下载地址http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/
fetchChromSizes hg19 > chrom_size.txt
### 提取对应的promoter区域
zcat gencode.v29lift37.annotation.gtf.gz | awk '$3=="gene" {print $0}' |  awk -v OFS="\t" '{if ($7=="+") {print $1, $4, $4+1} else {print $1, $5-1, $5}}' | bedtools slop -i - -g chrom_size.txt -b 5000 > promoter_5kb.bed

好啦，干货满满的介绍到此就结束了，希望大家有所收获。如果觉得该文章有帮助到你，欢迎转发，点好看，分享给更多的生信小伙伴。最后，再次感谢Ming Tang大神提供的资源。

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