如何計算每個基因的覆蓋度與深度

　如何計算每個基因的覆蓋度與深度，有多種方法可以完成。如下演示使用samtools depth命令方法

　　1. 數據下載

1.1 Fastq文件下載

　　從NCBI下載Illumina Hiseq X Ten平台的RNA-Seq數據SRR7751429信息如上圖所示。

1.1.1 使用wget命令（sra-toolkit工具下載太慢）下載

wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR775/SRR7751429/SRR7751429.sra

1.1.2 在SRA Toolkit工具頁面根據不同操作系統進行下載（例如，我的是編譯好的Centos 64位）

1.1.3 使用SRA toolkit工具將SRR7751429.sra數據轉成fastq格式

fastq-dump -split-3 SRR7751429.sra

1.2 基因組及注釋文件下載

　　人的參考基因組文件（版本GRCh38）下載

wget ftp://ftp.ensembl.org/pub/release-93/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz
# 解壓
gunzip Homo_sapiens.GRCh38.dna.toplevel.fa.gz

　　人的gtf注釋文件下載

wget ftp://ftp.ensembl.org/pub/release-93/gtf/homo_sapiens/Homo_sapiens.GRCh38.93.gtf.gz
# 解壓
gunzip Homo_sapiens.GRCh38.93.gtf.gz

　2. 生成bam文件

　　因為是演示，只需要生成bam文件，我這里就用bwa比對了，節約時間。

# 創建索引
bwa index /home/Ensembl/Animal/homo_sapiens/Homo_sapiens.GRCh38.dna.all.fa
# 比對
bwa mem -t 32 -M /home/Ensembl/Animal/homo_sapiens/Homo_sapiens.GRCh38.dna.all.fa SRR7751429_1.fastq SRR7751429_2.fastq -o SRR7751429.sam

　3. 基因CDS（編碼區）獲取

3.1 本地獲取基因cds信息

　　下載的Homo_sapiens.GRCh38.93.gtf文件包含有基因exon、cds、3'utr、5'utr等相關的物理位置信息，獲取基因CDS信息只需解析該文件就可以了。（有需要的話，后續跟新相關腳本）

3.2 使用ensembl獲取cds信息

　　

　　如上圖所示，以人BRCA2基因為例，搜到后點擊CCDS，出現該基因的物理位置信息。然后，將該信息復制粘貼，以如下圖所示格式儲存於文件BRCA2.bed中。

　4. 使用samtools工具進行統計

　　samtools工具是對SAM/BAM文件進行操作的軟件，其帶有多種統計相關的命令及SAM↔BAM格式轉換的命令。

4.1 SAM文件格式轉換為BAM文件格式

samtools view -@ 16 -bS SRR7751429.sam -o SRR7751429.bam

4.2 sort BAM文件，然后建立BAM文件索引

# sort BAM文件
samtools sort -@ 16 -o SRR7751429_sorted.bam SRR7751429.bam

# 索引BAM文件
samtools index -@ 16 SRR7751429_sorted.bam

4.3 使用depth命令計算bed文件區域中每個位點的深度

samtools depth  -b  BRCA2.bed SRR7751429_sorted.bam >BRCA2.bed.depth

　　一共得到3列以指標分隔符分隔的數據，第一列為染色體名稱，第二列為位點，第三列為覆蓋深度。

4.4 根據BED文件和深度文件來統計大於10×的區域占總CDS區域比例

# -*- coding: utf-8 -*-
from __future__ import division

import csv

# 定義cds文件名路徑
cdsfh = 'BRCA2.bed'

# 區域長度
cdslen = 0


with open(cdsfh, 'r') as f:
    cf = csv.reader(f, dialect='excel-tab')
    for row in cf:
        # 讀取每一行區域
        chrom, start, end = row
        length = int(end) - int(start) + 1
        # 迭代所有的cds區域長度，得到基因cds區域全長
        cdslen += length


# 定義深度文件名路勁
depthfh = 'BRCA2.bed.depth'

# 大於10X區域長度
gt10len = 0

with open(depthfh, 'r') as f:
    cf = csv.reader(f, dialect='excel-tab')
    for row in cf:
        # 讀取每一行區域
        chrom, pos, depth = row
        # 判斷覆蓋度是否大於10X,是的gt10len就自增1
        if int(depth) > 10: gt10len += 1


# 計算編碼區大於10X的區域占總編碼區的比例
percent = gt10len / cdslen * 100

# 輸出
print("%.2f%%" % percent)

　　上述腳本只能針對單個基因，若是多個基因，可結合shell循環實現。

參考資料

samtools

ensembl