bedtools神器 | gtf转bed | bed文件运算
我们生信技能书有一篇介绍bedtools的文章，可以在微信里搜着看下，非常有用。
http://bedtools.readthedocs.io/en/latest/
gtf转bed用Linux命令完全可以实现，因为gtf每一行比较规律，不像fasta和fastq。
cat gffcmp.combined.gtf | grep -v exon | cut -f1,4,5,9 | cut -f1 -d";" | awk '{print $1, $2, $3, $5}' | sed -e 's/ /\t/g' | sed -e 's/\"//g' > gffcmp.combined.bed
后面发现，只提取transcript的首尾是不对的，必须要提取exon信息：
cat gffcmp.combined.gtf | grep exon | cut -f1,4,5,9 | cut -f1 -d";" | awk '{print $1, $2, $3, $5}' | sed -e 's/ /\t/g' | sed -e 's/\"//g' > gffcmp.combined.exon.bed
awk轻松计算总的覆盖度：
cat cov.out.coding.exon.out | awk 'BEGIN{i=1;j=1}{i=i+$5;j=j+$6}END{print i,j}'
1. bed的 过滤/overlap 运算，就是从一个bed里过滤掉与另一个bed有交集的所有区域。
bedtools intersect -v -a gffcmp.combined.bed -b coding.bed -wa > test.bed
2. bed的覆盖度i计算，比如我有一些转录本，想知道它们在基因组上的覆盖度，怎么办，transcript之间是有overlap的
bedtools coverage -a Teatree_Assembly.fasta.bed -b gffcmp.combined.bed > cov.out
结果的每一行怎么解读，请看：链接
Sc0000000 1 3505831 977 1824630 3505830 0.5204559 Sc0000001 1 3488299 853 1929980 3488298 0.5532727 Sc0000002 1 2866215 896 1435119 2866214 0.5007020 Sc0000003 1 2774837 512 1285961 2774836 0.4634368 Sc0000004 1 2768176 402 720812 2768175 0.2603925
3. 2的另一种实现方式
bedtools genomecov -i Arabidopsis_thaliana.TAIR10.38.bed -g Arabidopsis_thaliana.TAIR10.dna_sm.bed > Arabidopsis_thaliana.cov
解读方法不一样：
chr1 0 980 1000 0.98 chr1 1 20 1000 0.02 chr2 1 500 500 1 genome 0 980 1500 0.653333 genome 1 520 1500 0.346667
chromosome (or entire genome) depth of coverage from features in input file number of bases on chromosome (or genome) with depth equal to column 2. size of chromosome (or entire genome) in base pairs fraction of bases on chromosome (or entire genome) with depth equal to column 2.
fai变bed
cat Arabidopsis_thaliana.TAIR10.dna_sm.fasta.fai | awk '{print $1, 1, $2}' | sed -e's/ /\t/g' > Arabidopsis_thaliana.TAIR10.dna_sm.bed
gff3转gtf
gffread Arabidopsis_thaliana.TAIR10.38.gff3 -T -o Arabidopsis_thaliana.TAIR10.38.gtf
gtf转bed
cat Arabidopsis_thaliana.TAIR10.38.gtf | cut -f1,4,5,9 | cut -f1 -d";" | awk '{print $1, $2, $3, $5}' | sed -e 's/ /\t/g' | sed -e 's/\"//g' | sed -e 's/transcript\://g' > Arabidopsis_thaliana.TAIR10.38.bed
coverage方法计算覆盖度
bedtools coverage -a Arabidopsis_thaliana.TAIR10.dna_sm.bed -b Arabidopsis_thaliana.TAIR10.38.coding.bed > coding.cov.c
最终统计
cat coding.cov.c | awk 'BEGIN{i=1;j=1}{i=i+$5;j=j+$6}END{print i,j}'
我比较了，这两种方法算出来的覆盖度是一摸一样的，但是如果真的只想算基因组覆盖度，genomecov肯定更快一些。
bedtools，真的神器，你值得拥有。
