周末閑著無(wú)聊在推特上閑逛拆内，刷著刷著發(fā)現(xiàn)了Ming Tang大神的一條推文。他將他日常經(jīng)常使用的單行命令整理到其Github中纤泵，然后分享給各位推友吮便。在點(diǎn)開(kāi)它的網(wǎng)站后笔呀，我猶如發(fā)現(xiàn)了一堆寶藏。里面很多單行命令行非常實(shí)用髓需，快速簡(jiǎn)單地解決了我日常所遇到的生物信息學(xué)數(shù)據(jù)處理問(wèn)題(沒(méi)有這些單行命令许师，我往往要通過(guò)自己寫(xiě)腳本或者借用其它工具來(lái)解決日常的生物信息學(xué)數(shù)據(jù)處理)。在取得Ming Tang大神的同意后授账，該推文將和大家整理分享這一系列實(shí)用的單行命令枯跑。
原文鏈接：
GitHub 地址： https://github.com/<wbr>crazyhottommy/bioinformatics-<wbr>one-liners
Fastq/Fasta文件的處理
統(tǒng)計(jì)fastq文件中的read的長(zhǎng)度與不同長(zhǎng)度的分布：
zcat file.fastq.gz | awk 'NR%4 == 2 {lengths[length($0)]++} END {for (l in lengths) {print l, lengths[l]}}'
fastq轉(zhuǎn)換為fasta格式:
zcat file.fastq.gz | paste - - - - | perl -ane 'print ">$F[0]\n$F[2]\n";' | gzip -c > file.fasta.gz
統(tǒng)計(jì)fastq文件中reads的數(shù)目:
cat file.fq | echo $((`wc -l`/4))
隨機(jī)抽取fastq文件的一部分，下面使用提取0.01%的reads作為例子(在不知道該方法之前白热，我一般使用Heng Li的seqtk sample來(lái)達(dá)到該目的):
cat file.fq | paste - - - - | awk 'BEGIN{srand(1234)}{if(rand() < 0.01) print $0}' | tr '\t' '\n' > out.fq
從fastq文件中提取含有某個(gè)堿基排列序列的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文件是同時(shí)包含了R1,R2的序列信息屋确，需要將它們分離出來(lái):
cat file.fq | paste - - - - - - - - | tee >(cut -f1-4 | tr '\t'
'\n' > out1.fq) | cut -f5-8 | tr '\t' '\n' > out2.fq
將一個(gè)含有多個(gè)染色體或者contigs的fasta文件纳击，按照染色體或者contigs分離成多個(gè)只含單一序列的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
反過(guò)來(lái)將單行的fasta文件，變成多行攻臀，下面例子對(duì)每60個(gè)堿基分一行:
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
統(tǒng)計(jì)mulitfasta中焕数，每個(gè)染色體/contigs的長(zhǎng)度:
awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen = seqlen +length($0)}END{print seqlen}' file.fa
對(duì)fasta文件中的header中添加對(duì)應(yīng)的信息:
### 目前你有的fasta文件
cat myfasta.txt
>Blap_contig79
MSTDVDAKTRSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>Bluc_contig23663
MSTNVDAKARSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>Blap_contig7988
MSTDVDAKTRSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
>Bluc_contig1223663
MSTNVDAKARSKERASIAAFYVGRNIFVTGGTGFLGKVLIEKLLRSCPDVGEIFILMRPKAGLSI
###需要添加的信息
cat my_info.txt
info1
info2
info3
info4
### 添加信息到對(duì)應(yīng)的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
簡(jiǎn)化fasta文件的表頭,比如將>7 dna:chromosome chromosome:GRCh37:7:1:159138663:1轉(zhuǎn)化為>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文件轉(zhuǎn)化成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文件轉(zhuǎn)化成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個(gè)CPU，按照染色體刨啸，并使用xargs去并行call variant能夠極大地提高其變異發(fā)現(xiàn)的效率:
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完變異后堡赔，也可以將不同染色體的變異合并起來(lái):
samtools view -H yourFile.bam | grep "\@SQ" | sed 's/^.*SN://g' | cut -f 1 | perl -ane 'system("cat tmp.$F[0].bcf >> yourFile.vcf");'
尋找當(dāng)前目錄下的bam文件，并使用5個(gè)CPUs將其復(fù)制到其它新的目錄中(這個(gè)命令在轉(zhuǎn)移文件時(shí)设联，非常實(shí)用):
find . -name "*bam" | xargs -P5 -I{} rsync -av {} dest_dir
根據(jù)vcf的表頭文件來(lái)sort VCF文件
cat my.vcf | awk '$0~"^#" { print $0; next } { print $0 | "sort -k1,1V -k2,2n" }'
將含有PASS的變異從VCF中提取出來(lái)(經(jīng)過(guò)一些工具過(guò)濾后善已，例如GATK灼捂，如果符合過(guò)濾的閥值，在VCF文件中會(huì)有一行PASS的字符):
cat my.vcf | awk -F '\t' '{if($0 ~ /\#/) print; else if($7 == "PASS") print}' > my_PASS.vcf
將bed文件按照不同的染色體分開(kāi)(同樣的道理可以應(yīng)用到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
處理其它文件/序列
將序列轉(zhuǎn)化其互補(bǔ)的序列换团，在設(shè)計(jì)primers很有用:
echo 'ATTGCTATGCTNNNT' | rev | tr 'ACTG' 'TGAC'
重命名一系列文件:
for file in *gz
do zcat $file > ${file/bed.gz/bed}
將看不到的字符顯示出來(lái)
### 方法一
cat my_file | sed -n 'l'
### 方法二
cat -A
復(fù)制大量的文件悉稠，使用rsync比copy速度更快
## 將所有文件從一個(gè)目錄復(fù)制到另一個(gè)目錄
rsync -av from_dir/ to_dir
##將所有文件從一個(gè)目錄復(fù)制到另一個(gè)目錄，但是已經(jīng)復(fù)制好的文件將會(huì)被跳過(guò):
rsync -avhP /from/dir /to/dir
顯示當(dāng)前目錄下所有的子目錄所占空間的大小
du -h --max-depth=1
查看當(dāng)前內(nèi)存的使用情況
free -mg
根據(jù)第一和第二列打印出特殊唯一的行:
awk '!a[$1,$2]++' input_file
### 或者
sort -u -k1,2 file
刪除空白行:
sed /^$/d
刪除最后一行
sed $d
使用!$,使用最后一個(gè)命令的字符
echo hello, world; echo !$
將最后執(zhí)行的命令保存為腳本:
echo "!!" > foo.sh
重復(fù)使用剛剛執(zhí)行命令的參數(shù):
!*
統(tǒng)計(jì)一個(gè)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}'
創(chuàng)建一個(gè)新的文件夾艘包，并且馬上cd進(jìn)去
mkdir blah && cd $_
切割所需要的行的猛，根據(jù)另一個(gè)文件提供的行名稱(chēng)(Thanks God,我一直想找現(xiàn)成的辦法去處理這個(gè)問(wèn)題,非常合適處理大型的文件，比R占更少的內(nèi)存):
#! /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}"
相同的問(wèn)題想虎，使用csvtk就簡(jiǎn)單很多了:
csvtk cut -t -f $(paste -s -d , list.txt) data.tsv
在每行開(kāi)始的地方增加或者刪除chr卦尊，有些軟件不支持chr字符，需要使用數(shù)字來(lái)表示染色體:
# 增加 chr
sed 's/^/chr/' my.bed
# 或者
awk 'BEGIN {OFS = "\t"} {$1="chr"$1; print}'
# 移除 chr
sed 's/^chr//' my.bed
檢查 tsv文件中行與列的數(shù)目的是否一致:
awk '{print NF}' test.tsv | sort -nu | head -n 1
只替代其中一列中的某個(gè)pattern
## 替代第五列
awk '{gsub(pattern,replace,$5)}1' in.file
## 使用強(qiáng)大的csvtk
csvtk replace -f 5 -p pattern -r replacement
反轉(zhuǎn)文件中的其中一列
###反轉(zhuǎn)第三列然后將其結(jié)果放到第五列中
awk -v OFS="\t" '{"echo "$3 "| rev" | getline $5}{print $0}'
###只反轉(zhuǎn)第二列
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可能超出了這一點(diǎn)邓线，這種情況下最好請(qǐng)使用bedtools slop根據(jù)基因組大小文件來(lái)避免這種情況。
## 使用fetchChromSizes提取染色體大小煌恢，下載地址http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/
fetchChromSizes hg19 > chrom_size.txt
### 提取對(duì)應(yīng)的promoter區(qū)域
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
好啦骇陈，干貨滿(mǎn)滿(mǎn)的介紹到此就結(jié)束了，希望大家有所收獲瑰抵。如果覺(jué)得該文章有幫助到你你雌，歡迎轉(zhuǎn)發(fā)，點(diǎn)好看二汛，分享給更多的生信小伙伴婿崭。最后，再次感謝Ming Tang大神提供的資源肴颊。
