2020-12-04生物信息單行腳本

生信單行腳本

生信息有用的單行腳本 (and some, more generally useful).

contents

Sources

Basic awk & sed

[返回頂部]

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

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

輸出第五列等于abc123的行:

awk '$5 == "abc123"' file.txt

輸出第五列不是abc123的行:

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

輸出第七列以字母a-f開頭的行:

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

輸出第七列不是以字母a-f開頭的行:

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

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

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

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

awk '$3>$5' file.txt

計算文件中第一列的累加值，輸出最后的結(jié)果:

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

計算第二列的平均值:

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

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

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

消除行開頭空和格制表符:

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

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

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

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

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

刪除空行:

sed '/^$/d' file.txt

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

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

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之間的所有行。（假設(shè)）染色體在第一列，位點在第三列（基于同樣的假設(shè)可以用來返回類似特定等位基因頻率的變異）

cat file.txt | awk '$1=="1"' | awk '$3>=1000000' | awk '$3<=2000000'

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:
基本序列統(tǒng)計。輸出總的reads數(shù)，不重復(fù)的reads總數(shù)，不重復(fù)reads百分比，最大冗余的序列及其頻度以及總占比百分?jǐn)?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}'

轉(zhuǎn)換.bam為.fastq:

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

Keep only top bit scores in blast hits (best bit score only):
只取blast采樣中的頂級位點的分?jǐn)?shù)（最高的位點分）

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

Keep only top bit scores in blast hits (5 less than the top):
只取blast采樣中的頂級位點的分?jǐn)?shù)（比頂級少于5的）

awk '{ if(!x[$1]++) {print $0; bitscore=($14-6)} 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; }'

轉(zhuǎn)化FASTQ文件為FASTA:

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

從第二行開始每四行取值（從FASTQ文件提取序列）。

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

輸出中剔除第一行：

awk 'NR>1' input.txt

輸出20-80行:

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

計算二，三行列的和并追加到每行后輸出

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

計算fastq文件平均reads的長度

awk 'NR%4==2{sum+=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,$2-1,$2,$4"/"$5,"+"}}'

sort, uniq, cut, etc.

[返回開頭]

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

cat -n file.txt

去重復(fù)行計數(shù)

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

找到兩文件都有的行（假設(shè)兩個文件都是無重復(fù)行，重定向執(zhí)行‘wd -l’計算同樣行的行數(shù)）

sort file1 file2 | uniq -d

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

# 用comm的方法
comm -12 file1 file2

對文件按照第九列數(shù)字順序排序（g按照常規(guī)數(shù)值，k列）

sort -gk9 file.txt

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

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

從文件中隨機(jī)取10行

shuf file.txt | head -n 10

輸出所有三個所可能的DNA序列

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

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的文件保存，這里假設(shè)/1 reads 在/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

Take a fasta file with a bunch of short scaffolds, e.g., labeled >Scaffold12345, remove them, and write a new fasta without them:

將一個fasta文件轉(zhuǎn)成一系列短的scaffolds。比如，標(biāo)簽 ">Scaffold12345"，然后移出他們，保存一個去掉他們的新文件：

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

Display hidden control characters:

顯示一個隱藏的控制字符：

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

find, xargs, and GNU parallel

[返回]

通過 https://www.gnu.org/software/parallel/. 載 GNU parallel

搜索文件夾及其子目錄中名稱為 .bam 文件（目錄也算）:

find . -name "*.bam"

刪除上面搜到的文件列表(不可逆的危險操作，謹(jǐn)慎使用！刪除之前請自習(xí)確認(rèn))

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

將所有.txt 文件修改為.bak(例如在對*.txt做操作之前用于文件備份)

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

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數(shù)據(jù)做Chastity過濾（grep 查詢 包含:N:，用（-A）選項第三列信息附加在匹配的包含一個序列質(zhì)量信息后，并保存為一個新的fasta文件）

find *fq | parallel "cat {} | grep -A 3 '^@.*[^:]*:N:[^:]*:' | grep -v '^\-\-$' > {}.filt.fq"

通過parallel并行運(yùn)行12個FASTQC任務(wù)

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

通過parallel給bam做索引，通過--dry-run打印測試這些命令，實際上并未做執(zhí)行。
find *.bam | parallel --dry-run 'samtools index {}'

seqtk

[back to top]

  • Seqtk項目托管地址https://github.com/lh3/seqtk。Seqtk是一個快捷輕量的處理FASTA和FASTQ格式基因序列的工具。他可以是先FASTA和FASTQ無縫處理和轉(zhuǎn)化，同時支持gzip格式的壓縮文件。

把FASTQ轉(zhuǎn)化為FASTA:

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

轉(zhuǎn)化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

折疊長FASTA/Q行，并且去除其注釋：

seqtk seq -Cl60 in.fa > out.fa

轉(zhuǎn)化多行FASTQ到四行FASTQ:

seqtk seq -l0 in.fq > out.fq

反轉(zhuǎn)FASTA/Q序列:

seqtk seq -r in.fq > out.fq

用序列文件中的名稱（比如name.1st）提取序列,一個虛列名一行:

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

利用序列文件中的”reg.bed“r信息提取地理信息的序列:

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

編碼‘reg.bed’地理信息到小寫

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

從兩個大的paired FASTQ文件提取10000個read pairs（記得用同樣的隨機(jī)種子保持 paire）

seqtk sample -s100 read1.fq 10000 > sub1.fq
seqtk sample -s100 read2.fq 10000 > sub2.fq

利用Phred公式從兩頭修剪低質(zhì)量bases:

seqtk trimfq in.fq > out.fq

從左端修剪5bp，從右端修剪10bp的。

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


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

GFF3 Annotations

[back to top]

輸出GFF3文件中標(biāo)注的所有的序列

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

檢測GFF3文件中標(biāo)注的所有性狀類型。

grep -v '^#' yourannots.gff3 | cut -s -f 3 | sort | uniq

檢測GFF3文件中標(biāo)注的基因數(shù)量。

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)'

FASTA頭列轉(zhuǎn)化為GFF格式（假設(shè)頭的長度，附加在”_length“ ,和Velvet assembled transcripts)）

grep '>' file.fasta | awk -F "_" 'BEGIN{i=1; print "##gff-version 3"}{ print $0"\t BLAT\tEXON\t1\t"$10"\t95\t+\t.\tgene_id="$0";transcript_id=Transcript_"i;i++ }' > file.gff

Other generally useful aliases for your .bashrc

有用的別名(.bashrc)

[back to top]

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

export PS1="\u@\h:\w\\$ "

避免反復(fù)敲諸如cd ../../..的命令（也可以用[autojump](https://github.com/joelthelion/autojump），讓你在飛速的轉(zhuǎn)換目錄

alias ..='cd ..'
alias ...='cd ../../'
alias ....='cd ../../../'
alias .....='cd ../../../../'
alias ......='cd ../../../../../'

向前和向后瀏覽

alias u='clear; cd ../; pwd; ls -lhGgo'
alias d='clear; cd -; ls -lhGgo'

覆蓋文件時候，先確認(rèn)

alias mv="mv -i"
alias cp="cp -i"  
alias rm="rm -i"

我最喜歡的”ls“別名

alias ls="ls -1p --color=auto"
alias l="ls -lhGgo"
alias ll="ls -lh"
alias la="ls -lhGgoA"
alias lt="ls -lhGgotr"
alias lS="ls -lhGgoSr"
alias l.="ls -lhGgod .*"
alias lhead="ls -lhGgo | head"
alias ltail="ls -lhGgo | tail"
alias lmore='ls -lhGgo | more'

對cut空格和逗號，分割文件

alias cuts="cut -d \" \""
alias cutc="cut -d \",\""

解壓縮tar包

alias tarup="tar -zcf"
alias tardown="tar -zxf"

或者可以用更普遍的‘extract’函數(shù)

# 源于ABSG(Advanced Bash Scripting Guide)中 Mendel Cooper的建議

extract () {
   if [ -f $1 ] ; then
       case $1 in
        *.tar.bz2)      tar xvjf $1 ;;
        *.tar.gz)       tar xvzf $1 ;;
        *.tar.xz)       tar Jxvf $1 ;;
        *.bz2)          bunzip2 $1 ;;
        *.rar)          unrar x $1 ;;
        *.gz)           gunzip $1 ;;
        *.tar)          tar xvf $1 ;;
        *.tbz2)         tar xvjf $1 ;;
        *.tgz)          tar xvzf $1 ;;
        *.zip)          unzip $1 ;;
        *.Z)            uncompress $1 ;;
        *.7z)           7z x $1 ;;
        *)              echo "don't know how to extract '$1'..." ;;
       esac
   else
       echo "'$1' is not a valid file!"
   fi
}

使用別名"mcd"創(chuàng)建一個目錄，并且cd到該目錄

function mcd { mkdir -p "$1" && cd "$1";}

跳轉(zhuǎn)到上級目錄，并且列出其內(nèi)容

alias u="cd ..;ls"

一個好看的grep

alias grep="grep --color=auto"

刷新你的.bashrc

alias refresh="source ~/.bashrc"

編輯你的.bashrc

alias eb="vi ~/.bashrc"

常用錯誤別稱

alias mf="mv -i"
alias mroe="more"
alias c='clear'

使用 pandoc轉(zhuǎn)化markdown文檔為PDF格式:

# 用法: mdpdf document.md document.md.pdf
alias mdpdf="pandoc -s -V geometry:margin=1in -V documentclass:article -V fontsize=12pt"

對當(dāng)前目錄搜索關(guān)鍵詞(ft "mytext" *.txt):

function ft { find . -name "$2" -exec grep -il "$1" {} \;; }

Etc

[返回]

重復(fù)運(yùn)行上一條命令:

sudo !!

列出最近最常用的命令行參數(shù)(通常是文件)

'ALT+.' or '<ESC> .'

敲出了部分命令，刪除這些輸入，查你忘記的明亮，拉回命令，繼續(xù)輸入(<CTRL+u>刪除光標(biāo)之前的輸入，<CTRL+y>恢復(fù)上個C-U刪除字符)

<CTRL+u> [...] <CTRL+y>

跳到一個目錄，執(zhí)行命令，然后返回當(dāng)前目錄(()的用法)

(cd /tmp && ls)

記時秒表 (輸入Enter or ctrl-d 停止):

time read

把上次執(zhí)行的命令生成一個腳本

echo "!!" > foo.sh

重用上次命令的所有參數(shù)

!*

列出或者刪除一個目錄中所有不匹配的特定后綴的文件（例如，列出所有不是壓縮的文件，刪除所有不以.foo和.bar后綴的文件）

ls !(*.gz)
rm !(*.foo|*.bar)

利用上次的命令，但是不需要他的的參數(shù)（重新輸入?yún)?shù)）:

!:- <new_last_argument>

激活一個快捷的編輯器，輸入，編輯長的，復(fù)雜，巧妙的命令:

fc

輸出一個特定的行（比如 42行）

sed -n 42p <file>

終結(jié)一個凍結(jié)的ssh session(會車換行，敲~鍵，在敲下.鍵）

[ENTER]~.

利用grep去除文件的空行，結(jié)果保存到新文件

grep . filename > newfilename

查找大文件（例如，大于500M的）

find . -type f -size +500M

利用截取列（例如，一個tab分割文件的第五個域）

cut -f5 --complement

查找包含特定字符的文件（-l 只輸出文件名, -i 忽略大小寫 -r 遍歷子目錄）

grep -lir "some text" *
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容