一些實用的生信單行命令

基本awk和sed命令

從文件中提取2,4,5列:

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

輸出第五列中包含abc123的行:

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

輸出第五列中不包含abc123的行:

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

輸出第七列中符合正則表達式規(guī)則(以a-f中任一字母開頭)的行:

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

輸出第七列中不符合正則表達式規(guī)則(以a-f中任一字母開頭)的行:

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

根據第二列去除重復行并輸出唯一行:

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

輸出第三列大于第五列的行:

awk '$3>$5' file.txt

對file.txt中的第一列求和:

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

計算第二列的均值:

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

將所有的foo替換為bar:

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

去除行首空格或制表符:

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

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

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

刪除空白行:

sed '/^$/d' file.txt

刪除包含EndOfUsefulData的行之后的所有內容:

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

維持原始順序同時刪除重復項:

awk '!visited[$0]++' file.txt

生信中實用的awk和sed命令

輸出1號染色體上處于1MB和2MB之間的行:

# 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)計fastq文件的一些基本信息,包括reads數量、唯一的reads數、唯一reads的比例、出現頻率最多的序列及頻數和所占比例:

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

從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

將vcf文件轉為bed文件:

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

根據reads名從fastq文件中提取指定reads:

zcat a.fastq.gz | awk 'BEGIN{RS="@";FS="\n"}; $1~/readsName/{print $2; exit}'

統(tǒng)計vcf中每行丟失的樣本數:

bcftools query -f '[%GT\t]\n' a.bcf |  awk '{miss=0};{for (x=1; x<=NF; x++) if ($x=="./.") {miss+=1}};{print miss}' > nmiss.count

sort, uniq, cut 命令

對文件中的每一行編號:

cat -n file.txt

統(tǒng)計文件中唯一行的數量:

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

查找兩個文件中相同的行:

sort file1 file2 | uniq -d

# Safer
sort -u file1 > a
sort -u file2 > b
sort a b | uniq -d

# Use comm
comm -12 file1 file2

按第九列進行數字排序:

sort -gk9 file.txt

查找第二列中出現最多的字符串:

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 兩個fastq文件:

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

find, xargs, 和 GNU parallel 命令

GNU parallel 下載地址:https://www.gnu.org/software/parallel/
在當前目錄中遞歸搜索bam文件:

find . -name "*.bam"

刪除所有的bam文件:

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

將所有的txt文件重命名成bak文件:

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 {}'

seqtk

Seqtk是一種快速處理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文件中所含區(qū)域的序列:

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

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

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

統(tǒng)計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)'

參考

?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容