统计bed文件下的reads数目和GC含量

最近发现,需要对bed(chr  start  end)文件进行处理,当你不知道有什么其他软件可以用的时候,不妨用bedtools来试试,看看它的各种用法,估计能满足。看名字,就大概知道干什么用了。

bedtools工具可用于广泛的基因组学分析任务,即基因组上的集合论。例如,bedtools允许人们在广泛使用的基因组文件格式(例如BAM,BED,GFF / GTF,VCF)中交叉,合并,计数,补充和混洗来自多个文件的基因组间隔。虽然每个单独的工具被设计为执行相对简单的任务(例如,交叉两个间隔文件),但是可以通过在UNIX命令行上组合多个bedtools操作来进行非常复杂的分析。

输入文件为:bam格式

例如:计算一个bed文件中的reads数和GC含量

bedtools map -a target.bed -b realigned.bam -c 10,10 -o count,concat | awk -v OFS="\t" '{n=length($5); gc=gsub("[gcGC]", "", $5); print $1,$2,$3,$4,gc/n}'>out.txt

其中reads数和bamdst处理的差不多,偶尔有几条的区别。

输入文件为:fasta格式

计算各碱基的数目

bedtools nuc [OPTIONS] -fi <fasta> -bed <bed/gff/vcf>


bedtools nuc -fi hg19.fa -bed target.bed |cut -f 1-4,5,13 >gc.txt results.txt

Output format: 
    The following information will be reported after each BED entry:
        1) %AT content
        2) %GC content
        3) Number of As observed
        4) Number of Cs observed
        5) Number of Gs observed
        6) Number of Ts observed
        7) Number of Ns observed
        8) Number of other bases observed
        9) The length of the explored sequence/interval.
        10) The seq. extracted from the FASTA file. (opt., if -seq is used)
        11) The number of times a user's pattern was observed.
            (opt., if -pattern is used.)
 

参考:https://www.plob.org/article/3748.html

        https://www.plob.org/article/10690.html

猜你喜欢

转载自blog.csdn.net/Cassiel60/article/details/89310811