搜索
查看: 2075|回复: 1

使用preseq计算文库复杂度以及估计加测量

[复制链接]

2

主题

2

帖子

65

积分

注册会员

Rank: 2

积分
65
发表于 2017-12-22 23:38:01 | 显示全部楼层 |阅读模式
本帖最后由 zxcippo 于 2017-12-23 11:00 编辑

在评估下机数据的时候,如果发现数据去重复之后无法达到目标覆盖度,那么就需要进一步加测。然而,有些文库复杂度很低,即使加测很多数据也无法得到更多的有效信息。那么如何评估文库复杂度,判断是否有加测的必要呢?

使用preseq软件可以实现根据现有测序数据评估已测序数据的复杂度,以及整个文库的复杂度。其中子命令  c_curve可以方便的计算现有测序数据中总测序量(total reads)与 有效数据量(distinct reads)的关系。流程如下:

下载安装 preseq
http://smithlabresearch.org/software/preseq/
安装教程中说是需要解压后 make all来安装,实际上直接解压后就可以执行。正常运行需要事先有samtools的环境变量。

使用preseq计算复杂度
[Bash shell] 纯文本查看 复制代码
preseq c_curve  -P -B  inbam > out.c_curve.txt

其中 -P代表是pairend数据,-B代表输入是bam(也可以使用bed作为输入文件,我没测试过)
得到的c_curve.txt 文件如下

如果有多个文件可以通过以下代码合并
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
# 2017-12-17
# zxcippo

$header = "gene";
undef %total;

foreach $i(0..@ARGV-1){
        $header .="\t$ARGV[$i]";
        open I,"$ARGV[$i]";
        while(<I>){
                chomp;
                @tmp = split;
                $total{$tmp[0]}[$i] = $tmp[1];
        }
}

print "$header\n";
foreach $key(keys %total){
        foreach $i (0..@ARGV-1){
                if(! $total{$key}[$i]){
                        $total{$key}[$i] ="NA";
                }
        }
        print $key,"\t",join("\t",@{$total{$key}}),"\n";
}



之后使用ggplot作图,参考代码如下

[Shell] 纯文本查看 复制代码
require(ggplot2)
require(reshape2)
require(stringr)
options(stringsAsFactors=F)
read.table("./merged.c_curve.txt",header=T)->ccurve
ccurve_melt<-melt(data=ccurve,id.vars="gene",variable.name = "sample",value.name="distinct_reads")
str_replace(ccurve_melt$sample,pattern=".c_curve.txt",replacement="")
colnames(ccurve_melt)[1]<-"total_reads"


ggplot(data=ccurve_melt) + geom_line(aes(x=total_reads,y=distinct_reads,group = sample,color=sample))
dev.off()


p<- ggplot(data=ccurve_melt) + geom_line(aes(x=total_reads,y=distinct_reads,group = sample)) + geom_abline(intercept=0,slope=1,lty=2) + ylim(0,12e07) + xlim(0,12e07) + theme(axis.text.x=element_text(angle=270)) + facet_wrap(~sample,nrow = 4)
ggsave(plot= p ,filename="estimate_complexity_single_mini_abline.pdf",device="pdf",units = "cm" , width= 25 ,height=25)


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:研究生第一篇论文常出现的9个问题
下一篇:SRA、SAM以及Fastq文件高速下载方法
回复

使用道具 举报

4

主题

56

帖子

597

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
597
发表于 2017-12-25 08:03:18 | 显示全部楼层
很不错,值得鼓励,看来又出现了一个比较不错的分享者
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-12 10:11 , Processed in 0.033002 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.