接著上一次內(nèi)容装蓬,繼續(xù)和大家分享常用實(shí)用的one-liners命令只嚣。
sort, uniq, cut, etc.
將文件中每一行的行號(hào)進(jìn)行標(biāo)識(shí):
cat -n file.txt
對(duì)文件中獨(dú)特的unique lines進(jìn)行計(jì)數(shù)(每周都會(huì)用到不下五次以上):
cat file.txt | sort -u | wc -l
對(duì)兩個(gè)單行的文件進(jìn)行sort穗椅,然后提取共有的元素(可以用來提取目標(biāo)基因ID):
sort -u file1 > a
sort -u file2 > b
comm -12 a b
對(duì)第幾行的數(shù)字進(jìn)行從小到大的排序:
sort -gk9 file.txt
隨機(jī)抽打碎文件的順序蔓倍。并抽取1000行 (個(gè)人用過來隨機(jī)抽取VCF文件中的變異位點(diǎn),來畫系統(tǒng)發(fā)育樹):
shuf file.txt | head -n 10
find, xargs, and GNU parallel
在目前的directory中遞歸地搜索以.bam 結(jié)尾的文件:
find . -name "*.bam"
刪除當(dāng)前目錄下所有以.bam結(jié)尾的文件(千萬要小心使用):
find . -name "*.bam" | xargs rm
將所有以.txt結(jié)尾的文件重新命名為將.bak 結(jié)尾的文件(一般用來backup一些文件敞咧,當(dāng)你想做一些大的處理或者改變的時(shí)候):
find . -name "*.txt" | sed "s/\.txt$//" | xargs -i echo mv {}.txt {}.bak | sh
這里會(huì)提到parallel漠魏,一個(gè)很強(qiáng)并行運(yùn)行文件的工具,下載鏈接:https://www.gnu.org/software/parallel/
同時(shí)并行運(yùn)行12個(gè)fastqc的jobs:
find *.fq | parallel -j 12 "fastqc {} --outdir ."
并行index你的bam file文件妄均,這里添加了一個(gè)(--dry-run
)的option柱锹,表示只會(huì)把命令打印出來,并不會(huì)真正的執(zhí)行丰包,可以用來給我檢查命令行有沒有按照我們的要求輸入對(duì):
find *.bam | parallel --dry-run 'samtools index {}'
seqtk
如果有了解過的同學(xué)禁熏,都會(huì)知道seqtk是一個(gè)很強(qiáng)大的Fasta/Fastq處理的一個(gè)工具,下載地址https://github.com/lh3/seqtk,或者直接通過conda install seqtk
去安裝邑彪。
將fastq格式的文件轉(zhuǎn)化成fasta格式
seqtk seq -a in.fq.gz > out.fa
提取name.list文件中含有的序列名稱的fa文件瞧毙,在name.list(也可以輸入bed文件)中一行有一個(gè)對(duì)應(yīng)的header名稱
seqtk subseq in.fa name.list > out.fa
從fastq文件中,隨機(jī)抽取10000 read:
seqtk sample -s100 read1.fq 10000 > sub1.fq
seqtk sample -s100 read2.fq 10000 > sub2.fq
將頭尾兩端低質(zhì)量的堿基切掉:
seqtk trimfq in.fq > out.fq
對(duì)Gff文件的處理
Gff文件也是一個(gè)我們?nèi)粘=?jīng)常遇到要處理的文件
找出你成功被注釋的序列名:
cut -s -f 1,9 yourannots.gff3 | grep $'\t' | cut -f 1 | sort | uniq
找出你gff文件中注釋的類型(exon/gene/mRNAd/CDS等)
grep -v '^#' yourannots.gff3 | cut -s -f 3 | sort | uniq
統(tǒng)計(jì)GFF文件中寄症,能注釋到gene的數(shù)目:
grep -c $'\tgene\t' yourannots.gff3
提取GFF文件中的gene ID:
grep $'\tgene\t' yourannots.gff3 | perl -ne '/ID=([^;]+)/ and printf("%s\n", $1)'
計(jì)算gff文件中宙彪,基因的長(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)'