Fork me on GitHub

WES数据分析流程(下)--可视化

字数统计: 7.1k字   |   阅读时长≈ 36分

[TOC]

转自Die Linke

MyNotes：

总体概括为：全外显子测序的标准分析流程通常包括以下步骤：

质量控制（Quality Control，QC）

对原始测序数据进行质量控制，包括去除低质量的序列片段、过滤掉低质量的碱基和低覆盖度的测序片段等。

读取比对（Read Alignment）

将过滤后的测序数据比对到参考基因组上，以确定每个片段的位置和方向。

变异检测（Variant Calling）

通过对比测序数据和参考基因组，识别出样本中存在的单核苷酸多态性（Single Nucleotide Polymorphism，SNP）和插入/删除多态性（Insertion/Deletion，InDel）等变异。

注释（Annotation）

对检测到的变异进行注释，包括位置信息、功能影响、已知疾病关联等。

过滤（Filtering）

根据预设的过滤标准，过滤掉可能为假阳性的变异，如常见的多态性、测序误差等。

数据解读（Data Interpretation）

将过滤后的变异进行进一步分析和解读，以确定与研究对象相关的变异。

功能验证（Functional Validation）

通过实验验证已识别的变异是否与疾病相关，例如验证变异对基因表达、蛋白质功能的影响等。

分析流程

  1. 准备工作:

  2. 数据上传:

从本地硬盘通过FileZilla软件或在云端通过scp协议将所需数据传输至进行分析的服务器.

  • FileZilla示例:

如下图所示，在窗口拖动数据文件到指定目录即可。

img

  • scp协议示例:

命令行输入以下命令进行文件传输。

1
2
3
4
5
6
scp -r local_folder remote_user_name@remote_ip:remote_folder
#-r表示递归传出路径下所有文件
#local_folder为传出的本地文件夹
#remote_user_name为接收文件的远程用户名
#remote_ip为接收文件的远程用户所在的ip地址
# remote_folder为远程用户下储存接收文件的文件夹
  1. 软件安装:
1
2
3
conda create -n WES_WGS_WGRS #新建一个名为WES_WGS_WGRS的conda环境
conda activate WES_WGS_WGRS#激活进入新建的WES_WGS_WGRS环境
conda install fastqc multiqc trimmomatic bwa samtools bcftools vcftools gffread #在该环境下安装fastqc等软件

img

  • 通过wget下载安装gtfToGenePred:
1
2
3
4
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64.v369/gtfToGenePred
mkdir -p /path/to/save/gtfToGenePred/bin #创建gtfToGenePred的储存路径并创建bin目录
cp gtfToGenePred /path/to/save/gtfToGenePred/bin #将gtfToGenePred移入储存路径下的bin目录
chmod +x /path/to/save/gtfToGenePred/bin/gtfToGenePred #赋予其读写执行全部权限
1
2
3
4
BiocManager::install('stringr')
BiocManager::install('circlize')
BiocManager::install('grid')
BiocManager::install('ComplexHeatmap')
  1. 参考基因组和基因组注释文件:

在对应的网站下载参考基因组和基因组注释文件, 如在NCBI的genome搜索ZM4能找到我要的物种对应的reference genome和annotation.

如图所示，点击红色方框内的genome和GFF会分别下载到.fasta格式的参考基因组序列和.gff格式的基因组注释。下方的基因组全长可以记录下来，在之后的分析中会使用到。

img

  1. 质量控制:

任何组学分析第一步都是对测序数据进行质量控制，细节可以参考之前的转录组教程 生科院的老张：RNAseq转录组差异表达分析教程 , 和转录组一样，使用md5检查数据完整性，fastqc进行质量检测，multiqc合并报告，trimmomatic进行质量过滤,参考代码如下：

  1. md5检查数据完整性：
1
2
3
4
md5sum *gz > md5.txt && md5sum -c md5.txt
# *gz 任意以gz结尾的文件
# > 运行结果保存至
# && 连接符，前一命令完成继续执行下一命令

2.2.2 fastqc质量控制与multiqc合并质控报告:

1
2
fastqc *gz && multiqc ./
# 对所有gz结尾文件进行质控并对当前目录下所有质控报告进行合并

2.2.3 trimmomatic质量过滤:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
trimmomatic 
PE #双端测序，单端测序为SE
-threads 4 #指定线程数为4
~/WES_WGS_WGRS/data/sample1/sample1_1.fq.gz #输入序列
~/WES_WGS_WGRS/data/sample1/sample1_2.fq.gz
~/WES_WGS_WGRS/data/sample1/sample1_paired_clean_1.fq.gz #输出配对序列和非配对序列
~/WES_WGS_WGRS/data/sample1/sample1_unpair_clean_1.fq.gz 
~/WES_WGS_WGRS/data/sample1/sample1_paired_clean_2.fq.gz
~/WES_WGS_WGRS/data/sample1/sample1_unpair_clean_2.fq.gz
ILLUMINACLIP:/data1/guest/yinlei/miniconda3/pkgs/trimmomatic-0.39-1/share/trimmomatic/adapters/TruSeq3-PE-2.fa:2:30:10:1:true 
#去除ILLUMINA接头，根据质控报告选择trimmomatic文件夹adapters路径下的接头文件
LEADING:3 #从reads开头切除质量低于阈值3的碱基
TRAILING:3 #从reads末尾切除质量低于阈值3的碱基
SLIDINGWINDOW:4:20 #从reads 5‘端开始进行长度为4的滑窗过滤，切除碱基质量低于阈值20的碱基
MINLEN:50 #丢弃剪切后长度低于阈值50的reads
TOPHRED33 #将reads的碱基质量值体系转为phred-33，若为phred-64则为TOPHRED64

由于一批测序数据中有多个测序文件，因此需要循环以上过程，参考代码如下图所示：

1
2
3
for dir in A B C D; #循环读取A,B,C,D四个目录
do XXXX; #执行XXXX
done #结束命令

img

trimmomatic质量过滤

  1. bwa序列比对:

将经过质量控制的序列与参考基因组进行比对，参考代码如下：

1
2
3
4
5
6
7
8
#建立比对索引
bwa index genome.fa

#将序列与基因组进行比对
bwa mem -t 16 genome.fa sample1_paired_clean_1.fq.gz sample1_paired_clean_2.fq.gz > sample1.sam
#参数如下：
# mem 比对算法
# -t 16 指定线程数为16

同样由于一批测序数据中有多个测序文件，因此需要循环以上过程，参考代码如下图所示：

img

bwa比对基因组

  1. samtools去除PCR重复后排序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#将比对结果sample1.sam转为.bam格式
samtools view -bS sample1.sam > sample1.bam
#参数如下:
# -bS 表示输出bam文件，输入为sam文件

#去除PCR重复
samtools rmdup -s sample1.bam sample1.rm.bam

#排序
samtools sort -@ 16 sample1.rm.bam > sample1.sorted.bam
#参数如下：
# -@ 16 指定线程数为16

#建立索引
samtools index sample1.sorted.bam

同样由于一批测序数据中有多个测序文件，因此需要循环以上过程，参考代码如下图所示：

img

samtools去除PCR重复，排序和建立索引

  1. samtools进行SNP calling：

得到基因组比对的去重复排序结果后，我们将比对信息中SNP的信息提取出来并整理成通用的突变格式VCF.

1
2
3
4
5
#SNP calling
samtools mpileup -gSDf genome.fa sample1.sorted.bam sample1.bcf

#将SNP鉴定结果sample1.bcf转为vcf格式
bcftools view -cvNg sample1.bcf sample1.vcf

同样由于一批测序数据中有多个测序文件，因此需要循环以上过程，参考代码如下图所示：

img

samtools进行SNP calling

img

bcftools进行vcf格式转换

  1. vcftools进行SNP质量过滤：

得到的VCF文件记录了我们分析得到的外显子(对应WES)或基因组(对应WGS或WGRS)上的全部SNP信息，但是这里边有一些是低质量的噪音信息，我们需要通过设置一系列的筛选标准来将噪音去除，称为SNP质量过滤或数据清洗。我一般选用以下三个标准：

  • –max-missing: 最大缺失值
  • –minDP: 最小测序深度
  • –minQ: 最小碱基质量

其中比较难以理解的是最大缺失值这个参数，vcftools官方给出的定义是这样的：

img

最大缺失值官方定义

我盯住的是0允许全部缺失而1代表不允许缺失，那就是说该值是介于0～1之间的浮点数，其大小代表着最少需要多少平行含有该SNP的比例。如0.3就代表至少30%的平行都有该SNP.

参考代码如下：

1
vcftools --vcf sample1.vcf --minDP 3 --max-missing 0.3 --minQ 30 --recode --recode-INFO-all --out sample1.filtered

同样由于一批测序数据中有多个测序文件，因此需要循环以上过程，参考代码如下图所示：

img

vcftools进行SNP质量过滤

  1. ANNOVAR进行SNP注释:

利用前面下载的基因组和注释文件对质量过滤后的SNP进行注释, 【将基因组和注释文件放在同一路径下】！

  1. gffread和gtfToGenePred进行注释文件的格式转换：
1
2
3
4
5
#我们下载的参考基因组注释文件是.gff格式的，需要转换为.gtf格式
gffread ZM4_annotation.gff -T -o ZM4_annotation.gtf
#将.gtf格式的注释文件转为ANNOVAR需要的输入文件
/path/to/bin/gtfToGenePred -genePredExt ZM4_annotation.gtf ZM4_refGene.txt
#【此时ZM4_refGene.txt应该和ZM4_annotation.gtf, ZM4_genome.fa是同一路径！！！】
  1. ANNOVAR进行注释：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
#建立注释索引
perl /path/to/annovar/retrieve_seq_from_fasta.pl --format refGene --seqfile ZM4_genome.fa ZM4_refGene.txt --outfile ZM4_refGeneMrna.fa
#【此时ZM4_refGeneMrna.fa应该和ZM4_refGene.txt, ZM4_annotation.gtf, ZM4_genome.fa是同一路径！！！】

#将SNP信息文件.vcf转为ANNOVAR需要的输入文件
perl /path/to/annovar/convert2annovar.pl -format vcf4old sample1.filtered.recode.vcf > sample1.annovar

#注释
perl /path/to/annovar/annotate_variation.pl -buildver ZM4 -outfile sample1.anno sample1.annovar /path/to/ZM4_refGeneMrna.fa and ZM4_refGene.txt/
#参数详解：
# -buildver ZM4 指定注释文件的版本是ZM4, 即注释文件XX_refGene.txt和XX_refGeneMrna.fa的前缀XX
# -outfile 输出注释文件
# 然后接输入文件
# 最后接注释文件的储存路径，即ZM4_refGeneMrna.fa, ZM4_refGene.txt, ZM4_annotation.gtf, ZM4_genome.fa的储存路径

SNP注释结果文件预览结果如下图所示：

img

由左至右各列信息为：

  • 所在行
  • 同义/非同义突变
  • 突变基因/外显子/核苷酸变化及位置/氨基酸变化及位置
  • 染色体
  • 突变起始位置
  • 突变终止位置
  • 原始序列
  • 突变序列
  • 同源/异源
  • FQ: Phred probability of all samples being the same
  • XXX
  • 测序质量

同样由于一批测序数据中有多个测序文件，因此需要循环注释过程，参考代码如下图所示：

img

ANNOVAR进行SNP注释

  1. 可视化补充信息—-samtools统计测序覆盖度和深度:

前方步骤中得到的去除PCR重复排序后的比对结果sample1.sorted.bam, 可通过samtools统计其中包含的测序覆盖度和深度信息用于可视化，参考代码如下：

1
samtools depth sample1.sorted.bam sample1.depth

得到的.depth文件预览结果如下：

  • 第一列为染色体名称
  • 第二列为染色体上的核苷酸位置
  • 第三列为该位置上的测序深度

img

samtools进行测序覆盖度和深度统计

可视化：

  1. 测序覆盖度和深度：

将上一步得到的.depth文件用于统计测序深度与覆盖度，在Python中的参考代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
import seaborn as sns
import matplotlib.pyplot as plt
######################################################PARAMETERS#######################################################
work_dir='/Users/ZYP/Downloads/WES_analysis_on_ZM4_RNAseq_data/3.raw/depth/' #输入和输出文件路径
files=['8P_MUT1','8P_OE1','8P_WT1','G5_MUT1','G5_OE1','G5_WT1','8P_MUT2','8P_OE2','8P_WT2','G5_MUT2','G5_OE2','G5_WT2'] #输入文件名前缀
start_position=254#from depth file 测序起始位置
end_position=2056015#from depth file 测序终止位置
distance=1000#distance for calculating average depth 用于统计平均深度和覆盖度的片段大小
length_genome=2056363#from genome 基因组总长
#######################################################################################################################
for name in files:
    list_pos = []
    list_depth = []
    sum_depths = 0
    distance_pos = 1
    for line in open(work_dir+name+'.depth'):
        sum_depths += float(line[:-1].split('\t')[2])
        if start_position >distance*distance_pos:
            list_pos.append(distance_pos)
            list_depth.append(0)
            distance_pos+=1
        if start_position==distance*distance_pos or start_position==length_genome:
            list_pos.append(distance_pos)
            if start_position==distance*distance_pos:
                list_depth.append(sum_depths/distance)
            else:
                list_depth.append(sum_depths/(length_genome%distance))
            sum_depths=0
            distance_pos += 1
        start_position+=1
    if start_position < length_genome:
        if start_position > distance*(distance_pos-1):
            list_pos.append(distance_pos)
            list_depth.append(sum_depths/(start_position-distance*(distance_pos-1)))
            distance_pos+=1
        list_pos.append(distance_pos)
        list_depth.append(0)
    print('data loading finished!')
    ax=sns.barplot(x=list_pos,y=list_depth)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['left'].set_color('gray')
    ax.spines['bottom'].set_color('gray')
    plt.savefig(work_dir+name+'_depth.pdf')
    plt.show()

运行结果如下图所示：

img

全基因组测序深度覆盖图

  1. SNP在基因组浏览器IGV的可视化:

前方得到的去重复排序后的比对结果文件.sorted.bam可以在基因组浏览器IGV (下载链接为 http://software.broadinstitute.org/software/igv/download ) 中打开，从而在基因组全局和单核苷酸局部进行个性化浏览和分析。

img

IGV可视化基因组

  1. SNV(单核苷酸变异)和InDel(插入缺失)分析：

包括以下部分：

  • SNV(单核苷酸变异，包括同义突变和非同义突变)
  • Nonsynonymous SNV(非同义突变单核苷酸变异)
  • Amino Acid Substitution(由非同义突变造成的氨基酸替换)

以ANNOVAR的SNP注释结果文件为输入，在Python中的参考代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
#######################################################PARAMETERS#######################################################
path='/Users/ZYP/Downloads/WES_analysis_on_ZM4_RNAseq_data/3.raw/annotation/' #输入和输出文件路径
files=['8P_MUT1','8P_OE1','8P_WT1','G5_MUT1','G5_OE1','G5_WT1','8P_MUT2','8P_OE2','8P_WT2','G5_MUT2','G5_OE2','G5_WT2'] #输入文件名前缀
########################################################################################################################
list_sample_SNV=[]
list_para_SNV=[]
list_SNV_name=[]
list_sample_SNV_nonsy=[]
list_para_SNV_nonsy=[]
list_SNV_nonsy_name=[]
list_sample_cAA=[]
list_para_cAA=[]
list_cAA_name=[]
save_file_cAA=open(path+'AA_substitution.txt','w')
save_file_cAA.write('Sample'+'\t'+'SNP type'+'\t'+'Gene'+'\t'+'REF NA'+'\t'+'ALT NA'+'\t'+'POS NA'+'\t'+'REF AA'+'\t'+'ALT AA'+'\t'+'POS AA'+'\n')
save_file_InDel=open(path+'InDel.txt','w')
save_file_InDel.write('Sample'+'\t'+'SNP type'+'\t'+'Hom/Het'+'\t'+'GeneFrom'+'\t'+'GeneIn'+'\t'+'Frame'+'\t'+'Effect'+'\n')

for file in files:
    sample=file[:-1]
    para=file[-1]
    dict_SNV = {'A>T|T>A':0,'A>C|T>G':0,'A>G|T>C':0,'C>G|G>C':0,'G>A|C>T':0,'G>T|C>A':0}
    dict_SNV_nonsy = {'A>T|T>A':0,'A>C|T>G':0,'A>G|T>C':0,'C>G|G>C':0,'G>A|C>T':0,'G>T|C>A':0}
    dict_cAA = {}
    for line in open(path+'SNP_annotation_'+file+'.txt'):
        kind=line.split('\t')[1]
        if kind in ['frameshift insertion','framshift deletion']:
            confusion=line.split('\t')[8]
            geneFrom_geneIn_frame_effect = line.split('\t')[2]
            geneFrom=geneFrom_geneIn_frame_effect.split(':')[0]
            geneIn=geneFrom_geneIn_frame_effect.split(':')[1]
            frame=geneFrom_geneIn_frame_effect.split(':')[-2]
            effect=geneFrom_geneIn_frame_effect.split(':')[-1]
            save_file_InDel.write(file+'\t'+kind+'\t'+confusion+'\t'+geneFrom+'\t'+geneIn+'\t'+frame+'\t'+effect+'\n')
        if kind in ['synonymous SNV','nonsynonymous SNV']:
            gene_cNA_cAA=line.split('\t')[2]
            gene=gene_cNA_cAA.split(':')[1]
            cNA=gene_cNA_cAA.split(':')[-2][2:]
            cAA=gene_cNA_cAA.split(':')[-1][2:-1]
            REF_NA=cNA[0]
            ALT_NA=cNA[-1]
            POS_NA=cNA[1:-1]
            REF_AA=cAA[0]
            ALT_AA=cAA[-1]
            POS_AA=cAA[1:-1]
            save_file_cAA.write(file+'\t'+kind+'\t'+gene+'\t'+REF_NA+'\t'+ALT_NA+'\t'+POS_NA+'\t'+REF_AA+'\t'+ALT_AA+'\t'+POS_AA+'\n')
            for key in dict_SNV.keys():
                list_changes=key.split('|')
                for change in list_changes:
                    if REF_NA==change[0] and ALT_NA==change[-1]:
                        dict_SNV[key]+=1
            if kind=='nonsynonymous SNV':
                name_cAA=REF_AA+'>'+ALT_AA
                for key in dict_SNV_nonsy.keys():
                    list_changes = key.split('|')
                    for change in list_changes:
                        if REF_NA==change[0] and ALT_NA==change[-1]:
                            dict_SNV_nonsy[key]+=1
                if name_cAA not in dict_cAA.keys():
                    dict_cAA[name_cAA]=1
                else:
                    dict_cAA[name_cAA]+=1
    for SNV in dict_SNV.keys():
        for time in range(0,dict_SNV[SNV]):
            list_sample_SNV.append(sample)
            list_para_SNV.append(para)
            list_SNV_name.append(SNV)
    for SNV_nonsy in dict_SNV_nonsy.keys():
        for time in range(0, dict_SNV[SNV_nonsy]):
            list_sample_SNV_nonsy.append(sample)
            list_para_SNV_nonsy.append(para)
            list_SNV_nonsy_name.append(SNV_nonsy)
    for cAA in dict_cAA.keys():
        for time in range(0, dict_cAA[cAA]):
            list_sample_cAA.append(sample)
            list_para_cAA.append(para)
            list_cAA_name.append(cAA)
data_SNV={'sample':list_sample_SNV,'parallel':list_para_SNV,'SNV':list_SNV_name}
data_SNV=pd.DataFrame(data_SNV)
ax1=sns.catplot(x='SNV', hue='parallel', col='sample', kind='count',data=data_SNV,palette='Set1')
plt.savefig(path+'SNV.pdf')
data_SNV_nonsy={'sample':list_sample_SNV_nonsy,'parallel':list_para_SNV_nonsy,'nonsynonymous SNV':list_SNV_nonsy_name}
data_SNV_nonsy=pd.DataFrame(data_SNV_nonsy)
ax2=sns.catplot(x='nonsynonymous SNV', hue='parallel', col='sample', kind='count',data=data_SNV_nonsy,palette='Set1')
plt.savefig(path+'Nonsynonymous_SNV.pdf')
data_cAA={'sample':list_sample_cAA,'parallel':list_para_cAA,'AA substitution':list_cAA_name}
data_cAA=pd.DataFrame(data_cAA)
ax3=sns.catplot(x='AA substitution', hue='parallel', col='sample', kind='count',data=data_cAA,palette='Set1')
plt.savefig(path+'AA_substitution.pdf')
save_file_cAA.close()

如图所示，运行结束得到可视化结果和记录SNV与InDel的两个文件SNV.txt和InDel.txt .

img

SNV and AA substitution

img

SNV and InDel

  1. 全基因组测序圈图可视化：

以前文中提到的测序深度文件.depth, ANNOVAR的SNP注释结果文件, 以及基因组注释文件.gff为输入，在R中的参考代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
library(stringr)    #方便处理字符串
library(circlize)    #绘制圈图
library(grid)    #设置画板
library(ComplexHeatmap)    #绘制图例

####################################PARAMETERS####################################
sample_name <-'G5_WT2'#测序样本名称
ref_name <-'ZM4_genome'#参考基因组名称
genome_size <-205636 #参考基因组长度
out_dir= '/Users/ZYP/Downloads/WES_analysis_on_ZM4_RNAseq_data/results/Analysis/' #输出路径
genome_gff <- '/Users/ZYP/Downloads/WES_analysis_on_ZM4_RNAseq_data/results/IGV visualization/ZM4_annotation.gff' #参考基因组gff注释文件
SNP_file <-'/Users/ZYP/Downloads/WES_analysis_on_ZM4_RNAseq_data/raw/annotation/SNP_annotation_G5_WT2.txt' #SNP注释文件
depth_base_stat <- '/Users/ZYP/Downloads/WES_analysis_on_ZM4_RNAseq_data/raw/depth/G5_WT2.depth' #测序深度覆盖度文件
seq_split <- 1000  #滑窗大小用于统计平均测序深度和SNP密度
##################################################################################

##测序覆盖度、平均深度统计
depth_base <- read.delim(depth_base_stat, stringsAsFactors = FALSE)
colnames(depth_base) <- c('gene','pos','depth')
depth <- depth_base$depth
depth_base <- depth_base[c(1,2)]
depth_base$end <- depth_base$pos
depth_base$depth <- depth

depth_exist <- subset(depth_base, depth != 0)
coverage <- round(100 * sum(depth_exist$seq_end - depth_exist$seq_start + 1) / genome_size, 2)
average_depth <- round(mean(depth_base$depth), 0)

depth_stat <- NULL
gene <- unique(depth_exist$gene)

for (gene_n in gene) {
  depth_subset <- subset(depth_exist, gene == gene_n)
  seq_end <- seq_split
  depth_num <- 0
  
  for (i in 1:nrow(depth_subset)) {
    if (depth_subset[i,'pos'] <= seq_end) depth_num <- depth_num + 1
    else {
      depth_stat <- rbind(depth_stat, c(gene_n, seq_end - seq_split + 1, seq_end,depth_num))
      
      seq_end <- seq_end + seq_split
      depth_num <- 0
      while (depth_subset[i,'pos'] > seq_end) {
        depth_stat <- rbind(depth_stat, c(gene_n, seq_end - seq_split + 1, seq_end, depth_num))
        seq_end <- seq_end + seq_split
      }
      depth_num <- depth_num + 1
    }
  }

  while (seq_end < depth_stat[gene_n,'seq_end']) {
    depth_stat <- rbind(depth_stat, c(gene_n, seq_end - seq_split + 1, seq_end, depth_num))
    seq_end <- seq_end + seq_split
    depth_num <- 0
  }
  depth_stat <- rbind(depth_stat, c(gene_n, seq_end - seq_split + 1, seq_stat[gene_n,'seq_end'], depth_num))
}

depth_stat <- data.frame(depth_stat, stringsAsFactors = FALSE)
names(depth_stat) <- c('gene', 'seq_start', 'seq_end', 'depth_num')
depth_stat$seq_start <- as.numeric(depth_stat$seq_start)
depth_stat$seq_end <- as.numeric(depth_stat$seq_end)
depth_stat$depth_num <- as.numeric(depth_stat$depth_num)
depth_stat

seq_stat <- NULL
for (gene in unique(depth_base$gene)) seq_stat <- rbind(seq_stat, c(gene, 0, 2056363))
seq_stat <- data.frame(seq_stat, stringsAsFactors = FALSE)
colnames(seq_stat) <- c('seq_ID', 'seq_start', 'seq_end')
rownames(seq_stat) <- seq_stat$seq_ID
seq_stat$seq_start <- as.integer(seq_stat$seq_start)
seq_stat$seq_end <- as.integer(seq_stat$seq_end)
seq_stat

##参考基因组基因信息，CDS & rRNA & tRNA
gene <- read.delim(genome_gff, header = FALSE, stringsAsFactors = FALSE, comment.char = '#')[c(1, 3, 4, 5, 7)]
gene <- subset(gene, V3 %in% c('CDS', 'rRNA', 'tRNA'))[c(1, 3, 4, 2, 5)]
names(gene) <- c('seq_ID', 'seq_start', 'seq_end', 'type', 'strand')

gene[which(gene$type == 'CDS'),'type'] <- 1
gene[which(gene$type == 'rRNA'),'type'] <- 2
gene[which(gene$type == 'tRNA'),'type'] <- 3
gene$type <- as.numeric(gene$type)
gene <- list(subset(gene, strand == '-')[-5], subset(gene, strand == '+')[-5])
gene

#读取 SNP 文件，统计 SNP 类型
snp <- read.delim(SNP_file, header = FALSE, colClasses = 'character', comment.char = '#')[c(4, 5, 7, 8)]
snp$V5 <- as.numeric(snp$V5)
snp$change <- str_c(snp$V7, snp$V8)

change <- which(snp$change == 'AT')
snp[change,'type1'] <- 'A>T|T>A'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'AG')
snp[change,'type1'] <- 'A>G|T>C'; snp[change,'type2'] <- 'ti'
change <- which(snp$change == 'AC')
snp[change,'type1'] <- 'A>C|T>G'; snp[change,'type2'] <- 'tv'

change <- which(snp$change == 'TA')
snp[change,'type1'] <- 'A>T|T>A'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'TG')
snp[change,'type1'] <- 'A>C|T>G'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'TC')
snp[change,'type1'] <- 'A>G|T>C'; snp[change,'type2'] <- 'ti'

change <- which(snp$change == 'GA')
snp[change,'type1'] <- 'G>A|C>T'; snp[change,'type2'] <- 'ti'
change <- which(snp$change == 'GT')
snp[change,'type1'] <- 'G>T|C>A'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'GC')
snp[change,'type1'] <- 'G>C|C>G'; snp[change,'type2'] <- 'tv'

change <- which(snp$change == 'CA')
snp[change,'type1'] <- 'G>T|C>A'; snp[change,'type2'] <- 'tv'
change <- which(snp$change == 'CT')
snp[change,'type1'] <- 'G>A|C>T'; snp[change,'type2'] <- 'ti'
change <- which(snp$change == 'CG')
snp[change,'type1'] <- 'G>C|C>G'; snp[change,'type2'] <- 'tv'

change_list <- c('A>T|T>A', 'A>G|T>C', 'A>C|T>G', 'G>A|C>T', 'G>T|C>A', 'G>C|C>G')
change_type <- c('synonymous SNV', 'nonsynonymous SNV')
anno_list <- list()
#length(anno_list) <- length(change_list) * length(change_type)
l = 0
anno <- read.delim(SNP_file, header = FALSE, colClasses = 'character', comment.char = '#')[c(4, 5, 5, 7, 12, 2)]
anno$V7 <-snp$type1
anno$V5 <- as.integer(anno$V5)
anno$V5.1 <- as.integer(anno$V5.1)
anno$V12 <- as.numeric(anno$V12)
colnames(anno)[1:6] <- c('seq_ID','seq_start', 'seq_end','change','ratio','type')
for (i in change_list) {
  for (j in change_type) {
    select <- anno[which(anno$change == i & anno$type == j), ]
    if (length(row.names(select))!=0) {
    l = l + 1
    anno_list[[l]] <- select
    }
  }
}

snp_ti <- length(which(snp$type2 == 'ti'))
snp_tv <- length(which(snp$type2 == 'tv'))

snp_at <- length(which(snp$type1 == 'A>T|T>A'))
snp_ag <- length(which(snp$type1 == 'A>G|T>C'))
snp_ac <- length(which(snp$type1 == 'A>C|T>G'))
snp_ga <- length(which(snp$type1 == 'G>A|C>T'))
snp_gt <- length(which(snp$type1 == 'G>T|C>A'))
snp_gc <- length(which(snp$type1 == 'G>C|C>G'))

#统计 SNP 密度
snp <- snp[c(1, 2, 5, 6, 7)]
colnames(snp)[1:2] <- c('seq_ID', 'seq_site')

snp_stat <- NULL
seq_ID <- unique(snp$seq_ID)

for (seq_ID_n in seq_ID) {
  snp_subset <- subset(snp, seq_ID == seq_ID_n)
  seq_end <- seq_split
  snp_num <- 0
  
  for (i in 1:nrow(snp_subset)) {
    if (snp_subset[i,'seq_site'] <= seq_end) snp_num <- snp_num + 1
    else {
      snp_stat <- rbind(snp_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, snp_num))
      
      seq_end <- seq_end + seq_split
      snp_num <- 0
      while (snp_subset[i,'seq_site'] > seq_end) {
        snp_stat <- rbind(snp_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, snp_num))
        seq_end <- seq_end + seq_split
      }
      snp_num <- snp_num + 1
    }
  }

  while (seq_end < seq_stat[seq_ID_n,'seq_end']) {
    snp_stat <- rbind(snp_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, snp_num))
    seq_end <- seq_end + seq_split
    snp_num <- 0
  }
  snp_stat <- rbind(snp_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_stat[seq_ID_n,'seq_end'], snp_num))
}

snp_stat <- data.frame(snp_stat, stringsAsFactors = FALSE)
names(snp_stat) <- c('seq_ID', 'seq_start', 'seq_end', 'snp_num')
snp_stat$seq_start <- as.numeric(snp_stat$seq_start)
snp_stat$seq_end <- as.numeric(snp_stat$seq_end)
snp_stat$snp_num <- as.numeric(snp_stat$snp_num)
snp_stat

#读取 SNP 文件，统计 InDel 长度
snp <- read.delim(SNP_file, header = FALSE, colClasses = 'character', comment.char = '#')[c(4, 5, 7, 8)]
snp$V5 <- as.numeric(snp$V5)
snp$change <- str_c(snp$V7, snp$V8)

indel <- read.delim(SNP_file, header = FALSE, colClasses = 'character', comment.char = '#')[c(4, 5, 7, 8)]
indel$V5 <- as.numeric(indel$V5)
indel$length <- str_length(indel[ ,4]) - str_length(indel[ ,3])
indel_insert <- length(which(indel$length > 0))
indel_delet <- length(which(indel$length < 0))

#统计 InDel 密度
indel <- indel[c(1, 2, 5)]
colnames(indel)[1:2] <- c('seq_ID', 'seq_site')

indel_stat <- NULL
seq_ID <- unique(indel$seq_ID)
for (seq_ID_n in seq_ID) {
  indel_subset <- subset(indel, seq_ID == seq_ID_n)
  seq_end <- seq_split
  indel_num <- 0
  
  for (i in 1:nrow(indel_subset)) {
    if (indel_subset[i,'seq_site'] <= seq_end) indel_num <- indel_num + 1
    else {
      indel_stat <- rbind(indel_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, indel_num))
      
      seq_end <- seq_end + seq_split
      indel_num <- 0
      while (indel_subset[i,'seq_site'] > seq_end) {
        indel_stat <- rbind(indel_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, indel_num))
        seq_end <- seq_end + seq_split
      }
      indel_num <- indel_num + 1
    }
  }
  
  while (seq_end < seq_stat[seq_ID_n,'seq_end']) {
    indel_stat <- rbind(indel_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_end, indel_num))
    seq_end <- seq_end + seq_split
    indel_num <- 0
  }
  indel_stat <- rbind(indel_stat, c(seq_ID_n, seq_end - seq_split + 1, seq_stat[seq_ID_n,'seq_end'], indel_num))
}

indel_stat <- data.frame(indel_stat, stringsAsFactors = FALSE)
names(indel_stat) <- c('seq_ID', 'seq_start', 'seq_end', 'indel_num')
indel_stat$seq_start <- as.numeric(indel_stat$seq_start)
indel_stat$seq_end <- as.numeric(indel_stat$seq_end)
indel_stat$indel_num <- as.numeric(indel_stat$indel_num)
indel_stat

pdf(str_c(out_dir, sample_name, '.circlize.pdf'), width = 40, height = 28)
circle_size = unit(1, 'snpc')
circos.par(gap.degree = 2, start.degree = 90)

#第一圈，染色体区域
circos.genomicInitialize(seq_stat, plotType = c('axis', 'labels'), major.by = 250000, track.height = 0.2)#, gap.degree=1

circos.genomicTrackPlotRegion(
  seq_stat, track.height = 0.05, stack = TRUE, bg.border = NA,
  panel.fun = function(region, value, ...) {
    circos.genomicRect(region, value, col = '#049a0b', border = NA, ...)
  } )

#第二圈，覆盖度 & 深度
circos.genomicTrack(
  depth_stat[1:4], track.height = 0.08,  ylim = c(0, (max(depth_stat$depth) + 1)), bg.col = '#EEEEEE6E', bg.border = NA,
  panel.fun = function(region,value, ...) {
    circos.genomicRect(region, value, ytop.column = 1, ybottom = 0, border = 'white', lwd = 0.2, col = 'red', ...)
    circos.lines(c(0, max(region)), c(average_depth, average_depth), col = 'red3', lwd=0.15, lty = 2)
    circos.yaxis( lwd = 0.1, tick.length = convert_x(0.15, 'mm'),labels.cex = 0.3)
  } )

#第三圈，CDS & rRNA & tRNA
color_assign <- colorRamp2(breaks = c(1, 2, 3), col = c('#00ADFF', 'orange', 'green2'))

circos.genomicTrackPlotRegion(
  gene, track.height = 0.16, stack = TRUE, bg.border = NA,
  panel.fun = function(region, value, ...) {
    circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
  } )

#第四圈，SNP 密度
value_max <- max(snp_stat$snp_num)
colorsChoice <- colorRampPalette(c('white', '#245B8E'))
color_assign <- colorRamp2(breaks = c(0:value_max), col = colorsChoice(value_max + 1))

circos.genomicTrackPlotRegion(
  snp_stat, track.height = 0.08, stack = TRUE, bg.border = NA,
  panel.fun = function(region, value, ...) {
    circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
  } )

#第五圈，InDel 密度
value_max <- max(indel_stat$indel_num)
colorsChoice <- colorRampPalette(c('white', '#7744A4'))
color_assign <- colorRamp2(breaks = c(0:value_max), col = colorsChoice(value_max + 1))

circos.genomicTrackPlotRegion(
  indel_stat, track.height = 0.08, stack = TRUE, bg.border = NA,
  panel.fun = function(region, value, ...) {
    circos.genomicRect(region, value, col = color_assign(value[[1]]), border = NA, ...)
  } )

#第六圈, SNV分布图
color_assign <- c('#BC80BD', '#FDB462', '#80B1D3', '#FB8072', '#8DD3C7', '#FFFFB3')
pch_assign <- c(16, 1)

circos.genomicTrackPlotRegion(
  anno_list, track.height = 0.08, bg.border = 'black', bg.lwd = 0.4,
  panel.fun = function(region, value, ...) {
  n = getI(...)
  circos.genomicPoints(region, value, cex = 1.0, pch = pch_assign[ifelse(n%%2, n%%2, 2)], col = color_assign[ifelse(n%%6, n%%6, 6)], ...)
  circos.yaxis(labels.cex = 0.3, lwd = 0.1, tick.length = convert_x(0.15, 'mm'))
  xlim = CELL_META$xlim
  ylim = CELL_META$ylim
  } )

#draw.sector(start.degree, end.degree, clock.wise = FALSE)#突出显示

#测序深度、覆盖度图例
depth_legend <- Legend(
  at = 1, labels = str_c(' Depth ( average: ', average_depth, ' X )'), labels_gp = gpar(fontsize = 8),
  title = str_c('Coverage: ', coverage , ' %'), title_gp = gpar(fontsize = 9),
  grid_height = unit(0.4, 'cm'), grid_width = unit(0.4, 'cm'), type = 'points', pch = NA, background = 'red')

#CDS & rRNA & tRNA 图例
gene_legend <- Legend(
  at = c(3, 2, 1), labels = c(' CDS', ' rRNA', ' tRNA'), labels_gp = gpar(fontsize = 8),
  title = 'CDS | rRNA | tRNA', title_gp = gpar(fontsize = 9), 
  grid_height = unit(0.4, 'cm'), grid_width = unit(0.4, 'cm'), type = 'points', pch = NA, background = c('#00ADFF', 'orange', 'green2'))

#SNP 密度图例
snp_legend <- Legend(
  at = round(seq(0, max(snp_stat$snp_num), length.out = 6), 0), labels_gp = gpar(fontsize = 8),
  col_fun = colorRamp2(round(seq(0, max(snp_stat$snp_num), length.out = 6), 0), colorRampPalette(c('white', '#245B8E'))(6)),
  title_position = 'topleft', title = 'SNP density', legend_height = unit(4, 'cm'), title_gp = gpar(fontsize = 9))

#SNV 分布图例
snv_legend <- Legend(
  at = c(1, 2, 3, 4, 5, 6, 7, 8),
  labels = c(' SNP: A>T|T>A', ' SNP: A>G|T>C', ' SNP: A>C|T>G', ' SNP: G>A|C>T', ' SNP: G>T|C>A', ' SNP: G>C|C>G', 'synonymous SNV', 'nonsynonymous SNV'),
  labels_gp = gpar(fontsize = 5), title = 'variant type\n', title_gp = gpar(fontsize = 6),
  grid_height = unit(0.4, 'cm'), grid_width = unit(0.4, 'cm'), type = 'points', background = NA,
  legend_gp = gpar(col = c('#BC80BD', '#FDB462', '#80B1D3', '#FB8072', '#8DD3C7', '#FFFFB3', 'gray', 'gray')),
  pch = c(20, 20, 20, 20, 20, 20, 16, 1) )

#左侧统计总览（涵括了上文大部分的统计概况，例如 SNP 替换类型统计等）
stat_legend <- Legend(
  at = 1, labels = '1', labels_gp = gpar(fontsize = 0), title_gp = gpar(fontsize = 9), 
  grid_height = unit(0, 'cm'), grid_width = unit(0, 'cm'), type = 'points', pch = NA, background = NA, 
  title = str_c('Sample: ', sample_name, '\nRefer species: ', ref_name, '\nRefer size: ', genome_size, ' bp\n\n\nTotal SNP: ', snp_ti + snp_tv, '\nTransitions: ', snp_ti, '\nTransversions: ', snp_tv, '\nTi/Tv: ', round(snp_ti / snp_tv, 2), '\nA>T|T>A: ', snp_at, '\nA>G|T>C: ', snp_ag, '\nA>C|T>G: ', snp_ac, '\nG>A|C>T: ', snp_ga, '\nG>T|C>A: ', snp_gt, '\nG>C|C>G: ', snp_gc, '\n\n\nTotal InDel: ', indel_insert + indel_delet, '\nInsert: ', indel_insert, '\nDelet: ', indel_delet))

#图例存放位置
y_coord <- 1.01
x_coord <- 0.275

pushViewport(viewport(x = x_coord+0.05, y = y_coord-0.208))
grid.draw(stat_legend)
upViewport()
pushViewport(viewport(x = x_coord-0.0005 , y = y_coord-0.165))
grid.draw(depth_legend)
upViewport()
pushViewport(viewport(x = x_coord-0.005 , y = y_coord-0.197))
grid.draw(gene_legend)
upViewport()
pushViewport(viewport(x = x_coord-0.0097 , y = y_coord-0.25))
grid.draw(snp_legend)
y_coord <- y_coord
upViewport()
pushViewport(viewport(x = x_coord-0.0097 , y = y_coord-0.3))
grid.draw(indel_legend)
upViewport()
pushViewport(viewport(x = 0.5, y = 0.5))
grid.draw(snv_legend)
upViewport()

circos.clear()
dev.off()#保存

运行结果如图所示，由外向内依次是：

  • 左上方图例和统计信息
  • 基因组
  • 染色体
  • 测序深度与覆盖度
  • 编码区和非编码区
  • SNV和InDel密度
  • SNV分布
  • SNV类型

img

全基因组测序圈图

img

全基因组测序圈图参数

扫一扫，分享到微信

tag:

    缺失模块。
    1、请确保node版本大于6.2
    2、在博客根目录（注意不是yilia-plus根目录）执行以下命令：
    npm i hexo-generator-json-content --save

    3、在根目录_config.yml里添加配置： 

      jsonContent:
    meta: false
    pages: false
    posts:
      title: true
      date: true
      path: true
      text: false
      raw: false
      content: false
      slug: false
      updated: false
      comments: false
      link: false
      permalink: false
      excerpt: false
      categories: false
      tags: true

主要涉及技术：
python R开发、开源爱好者、Linux

联系QQ:17769680069

志存高远

心怀善意