- 微信
- 微博
  
  分享文章到微博
- 复制链接
  
  复制链接到剪贴板

二代测序数据测序结果的比对和测序深度分析

benymorre 发表于 2019/01/05 15:37:16 2019/01/05

【摘要】二代测序数据一般是fastq 或 fastq.gz格式，分析流程大致如下： Fastq -> SAM -> BAM (排序，建索引) -> mpileup -> VCF -> VCF 的过滤和注释或者 Fastq -> SAM -> BAM (排序，建索引) -> VCF -> VCF 的过滤和注释以上两种分析流程，到BA...

二代测序数据一般是fastq 或 fastq.gz格式，分析流程大致如下：

Fastq -> SAM -> BAM (排序，建索引) -> mpileup -> VCF -> VCF 的过滤和注释

或者

Fastq -> SAM -> BAM (排序，建索引) -> VCF -> VCF 的过滤和注释

以上两种分析流程，到BAM这一步基本上是一致的，各个步骤用的软件会稍有差异。

---------------------------------------------------------------------

以下是用bwa进行比对产生bam的三步命令

bwa mem -R "@RG\tID:SR\tSM:NGS\tPL:Illunima" -p -t 8 /data/human/hg19.fasta test_1.fastq test_2.fastq > out.sam 2> align_withBWA.log
samtools view -@ 8 -bS out.sam | samtools sort - -@ 8 -T tmp_out.sorted -o out_sorted.bam
samtools index out_sorted.bam

Fastq 通过比对生成的文件就是SAM/BAM， BAM是二进制格式的SAM，处理起来速度更快些，占用的存储空间更小。

第三步samtools index 一定要做，因为很多处理BAM的软件都需要index文件，也就是后缀为bam.bai的文件

-----------------------------------------------------------------------

产生的BAM文件，这里用 GenomicAlignments 这个R包去分析

bam文件来自 samtools软件的测试bam（samtools-1.5/test/quickcheck/4.quickcheck.ok.bam"），我把文件名修改为out_sorted.bam

library(GenomicAlignments)
testAlign <- readGAlignments("./out_sorted.bam")

## 读取选定条件下的BAM

what <- c("flag", "cigar")

which <- GRanges('ref1', IRanges(8, 40))     ## 定义染色体编号和位点范围， ref1 如果是其它物种则可以是chr1 或者chrN等等
flag <- scanBamFlag(isMinusStrand = FALSE)
param <- ScanBamParam(which = which, what = what, flag=flag)

selectedBam <- readGAlignments("./out_sorted.bam", param=param)
cvg1 <- coverage(selectedBam)

## cvg1 中的测序深度可以直接查看， 在R中输入 cvg1 就可以看
## 如果要看在 ref1染色体  上的测序深度
cvg1$ref1

## 作图

plot(as.integer(cvg1$ref1),  pch=23, col='red', ylab='Coverage')

测序深度是指在基因组上的有关区域上，有多少个Reads覆盖到这些区域到各个位点，通常每个位点都有一个测序深度，最低是0.

结果如图：