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

在bioconductor上面看到一个R包 seqCNA ：

• PDF

• R Script

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

library(seqCNA) 
data(seqsumm_HCC1143) 
head(seqsumm_HCC1143)
dim(seqsumm_HCC1143)
tail(seqsumm_HCC1143)
## 200Kb windows to calculate the GC content and counts.
rco = readSeqsumm(tumour.data=seqsumm_HCC1143) 
rco = applyFilters(rco, trim.filter=1, mapq.filter=2)
rco = runSeqnorm(rco) 
rco = runGLAD(rco) 
plotCNProfile(rco) 
rco = applyThresholds(rco, seq(-0.8,4,by=0.8), 1)
plotCNProfile(rco)
summary(rco) 
head(rco@output)
writeCNProfile(rco,'./')

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

> head(seqsumm_HCC1143)
  chrom win.start reads.gc reads.mapq counts
1  chr1     0e+00    0.551      1.691   1199




    
2  chr1     2e+05    0.534      0.620    824
3  chr1     4e+05    0.457      0.831   8469
4  chr1     6e+05    0.545      6.479   1459
5  chr1     8e+05    0.720     31.619   1094
6  chr1     1e+06    0.766     38.205    976
> dim(seqsumm_HCC1143)
[1] 5314    5
> tail(seqsumm_HCC1143)
     chrom win.start reads.gc reads.mapq counts
5309  chr5 179800000    0.559     35.081   1946
5310  chr5 180000000    0.568     35.668   1970
5311  chr5 180200000    0.545     34.427   1790
5312  chr5 180400000    0.572     34.286   1569
5313  chr5 180600000    0.586     22.844   1591
5314  chr5 180800000    0.562      0.319    845

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

GENOME='/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta'
bam='ESCC13-T1_recal.bam'
samtools mpileup -f $GENOME $bam |\
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}' |\
sort 




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

得到的结果如下：

head GC_stat.10k.txt
chr1    1   7936    3885    582219
chr1    2   2123    934 88167
chr1    3   1969    169 28729
chr1    4   3556    593 48724
chr1    5   8828    2582    176627
chr1    6   8229    2290    117675
chr1    7   8794    438 156158
chr1    8   10000   723 211816
chr1    9   9077    2421    247285
chr1    10  9415    1661    371830

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

然后走seqCNA流程

library(seqCNA) 
a=read.table('wgs/GC_stat.10k.txt',fill = T)
a=na.omit(a)
a$GC = a[,4]/a[,3]
a$depth = a[,5]/a[,3]
a = a[a$depth<100,]
plot(a$GC,a$depth)
a=a[,c(1,2,6,5)]
colnames(a)=c('chrom','win.start','reads.gc','counts')
a$counts=floor(a$counts/150)
a$reads.mapq=30
library(seqCNA)  
head(a)
dim(a)
tail(a)
a$win.start=a$win.start*10000
a