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

單行命令俯逾,化繁為簡(jiǎn),重劍無(wú)鋒桌肴,大巧不工。接下來(lái)學(xué)習(xí)生信實(shí)用單行命令水醋。

基礎(chǔ)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開(kāi)頭的行:

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

輸出第7列中不以字母a-f開(kāi)頭的行:

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

計(jì)算第2列不重復(fù)的值保存在哈希arr中(一個(gè)值只保存一次):

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

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

awk '$3>$5' input.txt

計(jì)算文件第一列的累加值惶桐,輸出最后的結(jié)果:

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

計(jì)算第2列的平均值:

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

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

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

消除行首空格或制表符:

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

消除行結(jié)尾的空格和制表符:

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的片段之間的行信息潘懊,假設(shè)第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'

統(tǒng)計(jì)fastq文件的一些基本信息授舟,包括reads數(shù)、唯一的reads數(shù)岂却、唯一reads數(shù)的比例裙椭、出現(xiàn)頻率最多的序列及頻數(shù)和所占比例:(拿到rawdata來(lái)統(tǒng)計(jì)時(shí)很實(shí)用)

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文件轉(zhuǎn)fastq文件:

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

輸出blast結(jié)果中最高得分的結(jié)果:

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

將含多條序列fasta文件分割為多個(gè)fasta(每條一個(gè)文件):

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

輸出fasta文件中每條序列的序列名及其長(zhǎng)度:

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

將fastq文件轉(zhuǎn)為fasta文件:

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

從第2行開(kāi)始,每四行取值(從FASTQ文件提取序列):

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

輸出中剔出第1行:

awk 'NR>1' input.txt

輸出第20-80行:

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

計(jì)算第2扫尺、3行列的和,并追加到每行后輸出:

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

計(jì)算fastq文件評(píng)價(jià)read的長(zhǎng)度

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

轉(zhuǎn)化VSF文件為BED文件:

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

sort弊攘, uniq,cut

輸出帶行號(hào)的內(nèi)容:

cat -n file.txt

去除重復(fù)行并計(jì)數(shù):

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

找到兩文件都有的行(如果兩文件都是無(wú)重復(fù)行襟交,重定向執(zhí)行“wd -l”計(jì)算同樣行的行數(shù))

sort file1 file2 | uniq -d

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

# 用comm的方法
comm -12 file1 file2

文件按照第9列數(shù)字排序(g按照常規(guī)數(shù)值伤靠,k列):

sort -gk9 file.txt

找到第2列出現(xiàn)最多的字符:

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

從文件中隨機(jī)選取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兩個(gè)文件:(這里加上/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

將一個(gè)fasta文件轉(zhuǎn)成一系列短的scaffolds。比如:“>Scaffold12345”宴合,然后移除他們,保存一個(gè)去掉他們的新文件:

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

顯示一個(gè)隱藏的控制字符:

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等命令結(jié)合使用。
find:在指定目錄下查找蚤霞。

GNU parallel 下載地址:https://www.gnu.org/software/parallel/
搜素文件夾及其子目錄中名稱為.bam的文件(目錄也包括在內(nèi)):

find . -name "*.bam"

刪除所有的bam文件:(謹(jǐn)慎操作!U恪!)

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

將所有.txt文件重命名為.bak:(如在對(duì)*.txt做操作前用于文件備份)

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

同時(shí)運(yùn)行12個(gè)FASTQC文件:

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

將bam文件建索引(進(jìn)輸出命令滞乙,并不進(jìn)運(yùn)行程序):

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

seqktk

Seqtk專為FASTA和FASTQ而生,能快速處理FASTA或FASTQ格式序列序调,也可讀取gzip壓縮過(guò)的FASTA和FASTQ文件兔簇。下載地址:https://github.com/lh3/seqtk

將fastq轉(zhuǎn)為fasta格式:

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

將fastq文件中的質(zhì)量值低于20的序列屏蔽掉并轉(zhuǎn)為fasta格式:

# 將序列屏蔽為小寫
seqtk seq -aQ64 -q20 in.fq > out.fa
# 將序列屏蔽為N
seqtk seq -aQ64 -q20 -n N in.fq > out.fa

將fasta和fastq文件格式化為每行60個(gè)字符的多行序列并去除注釋信息:

seqtk seq -Cl60 in.fa > out.fa

將多行的fastq文件轉(zhuǎn)為4行的fastq文件:

seqtk seq -l0 in.fq > out.fq

生成fastq或fasta的反向互補(bǔ)序列:

seqtk seq -r in.fq > out.fq

根據(jù)name.lst(每行一個(gè)序列名)中的序列名提取序列:

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

提取reg.bed文件中所含區(qū)域的序列:

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

將reg.bed文件中所含區(qū)域的序列屏蔽為小寫:

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

使用Phred算法從兩端修剪低質(zhì)量的堿基:

seqtk trimfq in.fq > out.fq

從每條read的左端修剪5bp硬耍,從右端修剪10bp:

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

將合并的paired-end fastq文件拆分為-1和-2 兩個(gè)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

統(tǒng)計(jì)GFF3文件的基因數(shù)量:

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

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

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

統(tǒng)計(jì)GFF3文件每個(gè)基因的長(zhǎng)度

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

參考鏈接:

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末边酒,一起剝皮案震驚了整個(gè)濱河市,隨后出現(xiàn)的幾起案子坯认,更是在濱河造成了極大的恐慌,老刑警劉巖牛哺,帶你破解...
    沈念sama閱讀 212,542評(píng)論 6 493
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件劳吠,死亡現(xiàn)場(chǎng)離奇詭異,居然都是意外死亡痒玩,警方通過(guò)查閱死者的電腦和手機(jī),發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 90,596評(píng)論 3 385
  • 文/潘曉璐 我一進(jìn)店門燃观,熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái),“玉大人缆毁,你說(shuō)我怎么就攤上這事到涂〖箍颍” “怎么了践啄?”我有些...
    開(kāi)封第一講書人閱讀 158,021評(píng)論 0 348
  • 文/不壞的土叔 我叫張陵,是天一觀的道長(zhǎng)昭灵。 經(jīng)常有香客問(wèn)我,道長(zhǎng)烂完,這世上最難降的妖魔是什么诵棵? 我笑而不...
    開(kāi)封第一講書人閱讀 56,682評(píng)論 1 284
  • 正文 為了忘掉前任,我火速辦了婚禮履澳,結(jié)果婚禮上怀跛,老公的妹妹穿的比我還像新娘柄冲。我一直安慰自己吻谋,他們只是感情好羊初,可當(dāng)我...
    茶點(diǎn)故事閱讀 65,792評(píng)論 6 386
  • 文/花漫 我一把揭開(kāi)白布什湘。 她就那樣靜靜地躺著,像睡著了一般闽撤。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上贩据,一...
    開(kāi)封第一講書人閱讀 49,985評(píng)論 1 291
  • 那天闸餐,我揣著相機(jī)與錄音,去河邊找鬼舍沙。 笑死,一個(gè)胖子當(dāng)著我的面吹牛拂铡,可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播感帅,決...
    沈念sama閱讀 39,107評(píng)論 3 410
  • 文/蒼蘭香墨 我猛地睜開(kāi)眼,長(zhǎng)吁一口氣:“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼岖是!你這毒婦竟也來(lái)了?” 一聲冷哼從身側(cè)響起豺撑,我...
    開(kāi)封第一講書人閱讀 37,845評(píng)論 0 268
  • 序言:老撾萬(wàn)榮一對(duì)情侶失蹤硬梁,失蹤者是張志新(化名)和其女友劉穎,沒(méi)想到半個(gè)月后荧止,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體阶剑,經(jīng)...
    沈念sama閱讀 44,299評(píng)論 1 303
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡危号,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 36,612評(píng)論 2 327
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了猪半。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 38,747評(píng)論 1 341
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡磨确,死狀恐怖声邦,靈堂內(nèi)的尸體忽然破棺而出乏奥,到底是詐尸還是另有隱情亥曹,我是刑警寧澤,帶...
    沈念sama閱讀 34,441評(píng)論 4 333
  • 正文 年R本政府宣布媳瞪,位于F島的核電站,受9級(jí)特大地震影響蛇受,放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜龙巨,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 40,072評(píng)論 3 317
  • 文/蒙蒙 一旨别、第九天 我趴在偏房一處隱蔽的房頂上張望诗赌。 院中可真熱鬧秸弛,春花似錦、人聲如沸递览。這莊子的主人今日做“春日...
    開(kāi)封第一講書人閱讀 30,828評(píng)論 0 21
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽(yáng)镜雨。三九已至儿捧,卻和暖如春挑宠,著一層夾襖步出監(jiān)牢的瞬間,已是汗流浹背各淀。 一陣腳步聲響...
    開(kāi)封第一講書人閱讀 32,069評(píng)論 1 267
  • 我被黑心中介騙來(lái)泰國(guó)打工诡挂, 沒(méi)想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留碎浇,地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 46,545評(píng)論 2 362
  • 正文 我出身青樓璃俗,卻偏偏與公主長(zhǎng)得像奴璃,于是被迫代替她去往敵國(guó)和親。 傳聞我的和親對(duì)象是個(gè)殘疾皇子旧找,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 43,658評(píng)論 2 350

推薦閱讀更多精彩內(nèi)容