单个样本NGS数据如何做拷贝数变异分析呢

• PDF

• R Script

读其文档的时候发现，是可以针对单个样本进行拷贝数变异分析的，代码如下：

1. library(seqCNA)

2. data(seqsumm_HCC1143)

3. head(seqsumm_HCC1143)

4. dim(seqsumm_HCC1143)

5. tail(seqsumm_HCC1143)

6. ## 200Kb windows to calculate the GC content and counts.

7. rco = readSeqsumm(tumour.data=seqsumm_HCC1143)

8. rco = applyFilters(rco, trim.filter=1, mapq.filter=2)

9. rco = runSeqnorm(rco)

10. rco = runGLAD(rco)

11. plotCNProfile(rco)

12. rco = applyThresholds(rco, seq(-0.8,4,by=0.8), 1)

13. plotCNProfile(rco)

14. summary(rco)

15. head(rco@output)

16. writeCNProfile(rco,'./')

虽然其算法比较复杂，但是用法很简单，对input的数据进行了多步骤处理，而且其input数据本身是比较简单的，如下：

1. > head(seqsumm_HCC1143)

2.  chrom win.start reads.gc reads.mapq counts

3. 1  chr1     0e+00    0.551      1.691   1199

4. 2  chr1     2e+05    0.534      0.620    824

5. 3  chr1     4e+05    0.457      0.831   8469

6. 4  chr1     6e+05    0.545      6.479   1459

7. 5  chr1     8e+05    0.720     31.619   1094

8. 6  chr1     1e+06    0.766     38.205    976

9. > dim(seqsumm_HCC1143)

10. [1] 5314    5

11. > tail(seqsumm_HCC1143)

12.     chrom win.start reads.gc reads.mapq counts

13. 5309  chr5 179800000    0.559     35.081   1946

14. 5310  chr5 180000000    0.568     35.668   1970

15. 5311  chr5 180200000    0.545     34.427   1790

16. 5312  chr5 180400000    0.572     34.286   1569

17. 5313  chr5 180600000    0.586     22.844   1591

18. 5314  chr5 180800000    0.562      0.319    845

就是对单个样本的bam文件进行200kb的窗口进行滑动计算每个窗口的gc含量，该窗口区域覆盖的reads数量，还有比对的质量值，很容易写脚本进行计算。

1. GENOME='/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta'

2. bam='ESCC13-T1_recal.bam'

3. samtools mpileup -f $GENOME $bam |\

4. perl -alne '{$pos=int($F[1]/200000); $key="$F[0]\t$pos";$GC{$key}++ if $F[2]=~/[GC]/;$counts_sum{$key}+=$F[3];$number{$key}++;}END{print "$_\t$number{$_}\t$GC{$_}\t$counts_sum{$_}" foreach keys %number}' |\

5. sort -k1,1 -k 2,2n >T1.windows

得到的结果如下：

1. head GC_stat.10k.txt

2. chr1    1   7936    3885    582219

3. chr1    2   2123    934 88167

4. chr1    3   1969    169 28729

5. chr1    4   3556    593 48724

6. chr1    5   8828    2582    176627

7. chr1    6   8229    2290    117675

8. chr1    7   8794    438 156158

9. chr1    8   10000   723 211816

10. chr1    9   9077    2421    247285

11. chr1    10  9415    1661    371830

前面两行是窗口的坐标，第几号染色体的第几个窗口，后面3行是数据，分别是每个窗口的测到的碱基数，GC碱基数，测序总深度。

然后走seqCNA流程

1. library(seqCNA)

2. a=read.table('wgs/GC_stat.10k.txt',fill = T)

3. a=na.omit(a)

4. a$GC = a[,4]/a[,3]

5. a$depth = a[,5]/a[,3]

6. a = a[a$depth<100,]

7. plot(a$GC,a$depth)

8. a=a[,c(1,2,6,5)]

9. colnames(a)=c('chrom','win.start','reads.gc','counts')

10. a$counts=floor(a$counts/150)

11. a$reads.mapq=30

12. library(seqCNA)  

13. head(a)

14. dim(a)

15. tail(a)

16. a$win.start=a$win.start*10000

17. a=a[a$chrom %in% paste0('chr',c(1:22,'X','Y')),]

18. a=a[a$win.start>0,]

19. a=a[a$counts>0,]

20. a=a[a$reads.gc>0,]

21. a=a[,c(1:3,5,4)]

22. ## 200Kb windows to calculate the GC content and counts.

23. rco = readSeqsumm(tumour.data=a)

24. rco = applyFilters(rco, trim.filter=1, mapq.filter=2)

25. rco = runSeqnorm(rco)

26. rco = runGLAD(rco)

27. plotCNProfile(rco)

28. rco = applyThresholds(rco, seq(-0.8,4,by=0.8), 1)

29. plotCNProfile(rco)

30. summary(rco)

31. head(rco@output)

32. writeCNProfile(rco,'./')

分析结果汇总信息：

1. Basic information:

2.  SeqCNAInfo object with 306397 10Kbp-long windows.

3.  PEM information is not available.

4.  Paired normal is not available.

5.  Genome and build unknown (chromosomes chr1 to chrY).

6. Total filtered windows: 23717.

7. The profile is normalized and segmented.

8. Copy numbers called from 1 to 8.

9.  CN1:34924 windows.

10.  CN2:95226 windows.

11.  CN3:116829 windows.

12.  CN4:35047 windows.

13.  CN5:554 windows.

14.  CN6:100 windows.

15.  CN7:0 windows.

16.  CN8:0 windows.

单个样本的CNV结果如图

为什么要计算GC含量呢

这个是二代测序本身的技术限制，很容易探究到测序深度和GC含量是显著相关的，代码如下：

1. a=read.table('T1.windows')

2. a$GC = a[,4]/a[,3]

3. a$depth = a[,5]/a[,3]

4. a = a[a$depth<100,]

5. plot(a$GC,a$depth)

6. library(ggplot2)

7. # GET EQUATION AND R-SQUARED AS STRING

8. # SOURCE: http://goo.gl/K4yh

9. lm_eqn <- function(x,y){

10. m <- lm(y ~ x);

11. eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2,

12. list(a = format(coef(m)[1], digits = 2),

13. b = format(coef(m)[2], digits = 2),

14. r2 = format(summary(m)$r.squared, digits = 3)))

15. as.character(as.expression(eq));

16. }

17. p=ggplot(a,aes(GC,depth)) + geom_point() +

18. geom_smooth(method='lm',formula=y~x)+

19. geom_text(x = 0.5, y = 100, label = lm_eqn(a$GC , a$depth), parse = TRUE)

20. p=p+theme_set(theme_set(theme_bw(base_size=20)))

21. p=p+theme(text=element_text(face='bold'),

22. axis.text.x=element_text(angle=30,hjust=1,size =15),

23. plot.title = element_text(hjust = 0.5) ,

24. panel.grid = element_blank(),

25. #panel.border = element_blank()

26. )

27. print(p)

可以很明显看到GC含量和测序深度是高度相关的

GC含量和测序深度

WES的CNV探究-conifer软件使用