renluyao
/
quartet-wgs-germline-qc
Следить 1
В избранное 0
Форкнуть 0
Код Задачи 0 Pull Request'ы 0 Релизы 0 Вики Активность
Просмотр исходного кода

read me

master
LUYAO REN 4 лет назад
Родитель
2310dae4a7
Сommit
c25dae2952
2 измененных файлов: 83 добавлений и 193 удалений
Разделённый вид
Diff Options
Show Stats Download Patch File Download Diff File
  1. +66
    -178
      README.md
  2. +17
    -15
      codescripts/extract_tables.py

+ 66
- 178
README.md Просмотреть файл

@@ -160,126 +160,30 @@ Date是指数据获得日期，格式为20200710

## App输出文件

本计算会产生大量的中间结果，这里说明最后整合好的结果文件，上述的三种情况得到的结果文件如下。
本计算会产生大量的中间结果，这里说明最后整合好的结果文件。两个tasks输出最终的结果：

###1. 没有测LCL5和LCL6，或者没有同时测LCL5和LCL6
#### 1. extract_tables.wdl

####1.1 原始数据质量控制
原始结果质控 pre_alignment.txt

####1.2 比对后数据质量控制
比对结果指控 post_alignment.txt

##### 输出目录
突变检出指控 variants.calling.qc.txt

extract_multiqc
如果用户输入4个一组完整的家系样本则可以得到每个家庭单位的precision和recall的平均值，用于（1）计算综合指标和（2）报告第一页的展示：

##### 输出结果文件
precision.recall.txt

qualimap.final.result.txt
####2. quartet_mendelian.wdl

####1.3 突变检出数据质量控制
基于Quartet家系的质控 ${project}.mendelian.txt

#####与标准集进行比较
#### 3. D5_D6.WDL

######输出目录

extract_multiqc

######输出结果文件

benchmark.final.result.txt

###2. 包含LCL5和LCL6同卵双胞胎的数据，但是父母的数据不全

#### 2.1 原始数据质量控制

##### 输出目录

extract_multiqc

##### 输出结果文件

fastqc.final.result.txt

fastqscreen.final.result.txt

#### 2.2 比对后数据质量控制

##### 输出目录

extract_multiqc

##### 输出结果文件

qualimap.final.result.txt

#### 2.3 突变检出数据质量控制

##### 2.3.1 与标准集进行比较

#####输出目录

extract_multiqc

#####输出结果文件

benchmark.final.result.txt

**2.3.2** **姐妹一致性**

##### 输出目录

D5_D6

##### 输出结果文件
如果用户没有完整输入一组家庭，但有同时有D5和D6的信息，我们可以计算同卵双胞胎检测出的突变一致性，但是这部分输出暂不整合至报告中。

${project}.sister.txt

### 3. 四个quartet样本都测了

#### 3.1 原始数据质量控制

##### 输出目录

extract_multiqc

##### 输出结果文件

fastqc.final.result.txt

fastqscreen.final.result.txt

#### 3.2 比对后数据质量控制

##### 输出目录

extract_multiqc

##### 输出结果文件

qualimap.final.result.txt

#### 3.3 突变检出数据质量控制

##### 3.3.1 与标准集进行比较

###### 输出目录

extract_multiqc

###### 输出结果文件

benchmark.final.result.txt

##### 3.3.2 通过Quartet家系设计

###### 输出目录

quartet_mendelian

###### 输出结果文件

${project}.mendelian.txt

## 结果展示与解读

####1. 原始数据质量控制
@@ -288,83 +192,67 @@ ${project}.mendelian.txt

FastQC和FastqScreen是两个常用的原始数据质量控制软件

#####1.1 fastqc.final.result.txt

| 列名 | 说明 |
| -------------------------------------------------- | ------------------------------------------------------------ |
| Sample | 样本名，R1结尾为read1，R2结尾为read2 |
| FastQC_mqc-generalstats-fastqc-percent_duplicates | % Duplicate reads |
| FastQC_mqc-generalstats-fastqc-percent_gc | Average % GC content |
| FastQC_mqc-generalstats-fastqc-avg_sequence_length | Reads长度 |
| FastQC_mqc-generalstats-fastqc-percent_fails | 12个Fastqc模块中失败的比例 |
| FastQC_mqc-generalstats-fastqc-total_sequences | Total sequences |
| per_base_sequence_quality | pass/warn/fail；Quality就是Fred值，-10*log10(p)，p为测错的概率。所以一条reads某位置出错概率为0.01时，其quality就是20。横轴代表位置，纵轴quality。红色表示中位数，黄色是25%-75%区间，触须是10%-90%区间，蓝线是平均数。若任一位置的下四分位数低于10或中位数低于25，报"WARN"；若任一位置的下四分位数低于5或中位数低于20，报"FAIL"。 |
| per_tile_sequence_quality | pass/warn/fail；每个tile的测序情况，横轴是测序序列第1个碱基到第101个碱基，纵轴是tail的Index编号，这个图主要是为了防止，在测序过程中，某些tail受到不可控因素的影响而出现测序质量偏低，蓝色代表测序质量很高，暖色代表测序质量不高，如果某些tail出现暖色，可以在后续分析中把该tail测序的结果全部都去除。 |
| per_sequence_quality_scores | pass/warn/fail；每条reads的quality的均值的分布，横轴为quality，纵轴是reads数目。当出现上图的情况时，我们就会知道有一部分reads具有比较差的质量。当峰值小于27（错误率0.2%）时报"WARN"，当峰值小于20（错误率1%）时报"FAIL"。 |
| per_base_sequence_content | pass/warn/fail；对所有reads的每一个位置，统计ATCG四种碱基（正常情况）的分布：横轴为位置，纵轴为百分比。 正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。当所有位置的碱基比例一致的表现出bias时，即四条线平行但分开，往往代表文库有bias (建库过程或本身特点)，或者是测序中的系统误差。当任一位置的A/T比例与G/C比例相差超过10%，报"WARN"；当任一位置的A/T比例与G/C比例相差超过20%，报"FAIL"。 |
| per_sequence_gc_content | pass/warn/fail；统计reads的平均GC含量的分布。红线是实际情况，蓝线是理论分布（正态分布，均值不一定在50%，而是由平均GC含量推断的）。 曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差（overrepresented reads）。形状接近正态但偏离理论分布的情况提示我们可能有系统偏差。偏离理论分布的reads超过15%时，报"WARN"；偏离理论分布的reads超过30%时，报"FAIL"。 |
| per_base_n_content | pass/warn/fail；当测序仪器不能辨别某条reads的某个位置到底是什么碱基时，就会产生“N”。对所有reads的每个位置，统计N的比率：正常情况下N的比例是很小的，所以图上常常看到一条直线，但放大Y轴之后会发现还是有N的存在，这不算问题。当Y轴在0%-100%的范围内也能看到“鼓包”时，说明测序系统出了问题。当任意位置的N的比例超过5%，报"WARN"；当任意位置的N的比例超过20%，报"FAIL"。 |
| sequence_length_distribution | pass/warn/fail；reads长度的分布。当reads长度不一致时报"WARN"；当有长度为0的read时报“FAIL”。 |
| sequence_duplication_levels | pass/warn/fail；统计序列完全一样的reads的频率。测序深度越高，越容易产生一定程度的duplication，这是正常的现象，但如果duplication的程度很高，就提示我们可能有bias的存在（如建库过程中的PCR duplication）。横坐标是duplication的次数，纵坐标是duplicated reads的数目，以unique reads的总数作为100%。 如果原始数据很大（事实往往如此），做这样的统计将非常慢，所以fastqc中用fq数据的前200,000条reads统计其在全部数据中的重复情况。重复数目大于等于10的reads被合并统计，这也是为什么我们看到上图的最右侧略有上扬。大于75bp的reads只取50bp（不知道怎么选的）进行比较。但由于reads越长越不容易完全相同（由测序错误导致），所以其重复程度仍有可能被低估。当非unique的reads占总数的比例大于20%时，报"WARN"；当非unique的reads占总数的比例大于50%时，报"FAIL“。 |
| overrepresented_sequences | pass/warn/fail；如果有某个序列大量出现，就叫做over-represented。fastqc的标准是占全部reads的0.1%以上。和上面的duplicate analysis一样，为了计算方便，只取了fq数据的前200,000条reads进行统计，所以有可能over-represented reads不在里面。而且大于75bp的reads也是只取50bp。如果命令行中加入了-c contaminant file，出现的over-represented sequence会从contaminant_file里面找匹配的hit（至少20bp且最多一个mismatch），可以给我们一些线索。当发现超过总reads数0.1%的reads时报”WARN“，当发现超过总reads数1%的reads时报”FAIL“。 |
| adapter_content | pass/warn/fail；序列中两端adapter的情况；如果有adapter序列没有去除干净的情况，在后续分析的时候需要先使用cutadapt软件进行去接头，也可以用 trimmomatic来去除接头。 |
| kmer_content | pass/warn/fail；如果某k个bp的短序列在reads中大量出现，其频率高于统计期望的话，fastqc将其记为over-represented k-mer。默认的k = 5，可以用-k --kmers选项来调节，范围是2-10。出现频率总体上3倍于期望或是在某位置上5倍于期望的k-mer被认为是over-represented。fastqc除了列出所有over-represented k-mers，还会把前6个的per base distribution画出来。当有出现频率总体上3倍于期望或是在某位置上5倍于期望的k-mer时，报”WARN“；当有出现频率在某位置上10倍于期望的k-mer时报"FAIL"。 |

#####1.2 fastqscreen.final.result.txt

| 列名 | 说明 |
| ------------------ | ------------------------------------------------------------ |
| Sample | 样本名，R1结尾为read1，R2结尾为read2 |
| Human percentage | 比对到人类基因组的比例 |
| ERCC percentage | 比对到External RNA Controls Consortium基因组的比例 |
| EColi percentage | 比对到大肠杆菌基因组的比例 |
| Adapter percentage | 比对到接头序列的比例 |
| Vector percentage | 比对到载体基因组的比例 |
| rRNA percentage | 比对到rRNA序列的比例 |
| Virus percentage | 比对到病毒基因组的比例 |
| Yeast percentage | 比对到酵母基因组的比例 |
| Mitoch percentage | 比对到线粒体序列的比例 |
| Phix percentage | 比对到Phix基因组的比例。PhiX对照品v3是一款可靠、连接接头的文库，适合用作Illumina测序运行的对照品。 |
| No hits percentage | 没有比对到以上基因组的比例 |
本APP的输出文件是 **pre_alignment.txt**

| 列名 | 说明 |
| ------------------------- | ------------------------------------ |
| Sample | 样本名，R1结尾为read1，R2结尾为read2 |
| %Dup | % Duplicate reads |
| %GC | Average % GC content |
| Total Sequences (million) | Total sequences |
| %Human | 比对到人类基因组的比例 |
| %EColi | 比对到大肠杆菌基因组的比例 |
| %Adapter | 比对到接头序列的比例 |
| %Vector | 比对到载体基因组的比例 |
| %rRNA | 比对到rRNA序列的比例 |
| %Virus | 比对到病毒基因组的比例 |
| %Yeast | 比对到酵母基因组的比例 |
| %Mitoch | 比对到线粒体序列的比例 |
| %No hits | 没有比对到以上基因组的比例 |

#### 2. 比对后数据质量控制

##### qualimap.final.result.txt

| 列名 | 说明 |
| ----------------------------------------------------- | --------------------------------------------- |
| Sample | 样本名 |
| QualiMap_mqc-generalstats-qualimap-avg_gc | Mean GC content |
| QualiMap_mqc-generalstats-qualimap-median_insert_size | Median insert size（bp） |
| QualiMap_mqc-generalstats-qualimap-1_x_pc | Fraction of genome with at least 1x coverage |
| QualiMap_mqc-generalstats-qualimap-5_x_pc | Fraction of genome with at least 5x coverage |
| QualiMap_mqc-generalstats-qualimap-10_x_pc | Fraction of genome with at least 10x coverage |
| QualiMap_mqc-generalstats-qualimap-30_x_pc | Fraction of genome with at least 30x coverage |
| QualiMap_mqc-generalstats-qualimap-50_x_pc | Fraction of genome with at least 50x coverage |
| QualiMap_mqc-generalstats-qualimap-median_coverage | Median coverage |
| QualiMap_mqc-generalstats-qualimap-percentage_aligned | % mapped reads |
| QualiMap_mqc-generalstats-qualimap-mapped_reads | Number of mapped reads |
| QualiMap_mqc-generalstats-qualimap-total_reads | Total reads |
| QualiMap_mqc-generalstats-qualimap-general_error_rate | Mapping error rate |
本APP的输出文件是 **post_alignment.txt**

| 列名 | 说明 |
| --------------------- | --------------------------------------------- |
| Sample | 样本名 |
| %Mapping | % mapped reads |
| %Mismatch Rate | Mapping error rate |
| Mendelian Insert Size | Median insert size（bp） |
| %Q20 | % bases >Q20 |
| %Q30 | % bases >Q30 |
| Median Coverage | Median deduped coverage |
| PCT_1X | Fraction of genome with at least 1x coverage |
| PCT_5X | Fraction of genome with at least 5x coverage |
| PCT_10X | Fraction of genome with at least 10x coverage |
| PCT_30X | Fraction of genome with at least 30x coverage |

####3. 突变检出数据质量控制

##### 3.1 与标准集对比 benchmark.final.result.txt

| 列名 | 说明 |
| --------- | ------ |
| Sample | 样本名 |
| Precision | 查准率 |
| Recall | 查全率 |

##### 3.2 姐妹一致性 ${project}.sister.txt （样本没有没有测全，但是同时测了LCL5和LCL6）

| 列名 | 说明 |
| --------------------- | ------------------------------------------------------------ |
| Family | 家庭名字，我们目前的设计是4个Quartet样本，每个三个技术重复，family_1是指rep1的4个样本组成的家庭单位，以此类推。 |
| Reproducibility_D5_D6 | Quartet-D5和Quartet-D6的一致性 |

#####3.3 Quartet家系关系评估 ${project}.mendelian.txt
本APP的输出文件是 **variants.calling.qc.txt**

| 列名 | 说明 |
| --------------- | ------------------------------ |
| Sample | 样本名 |
| SNV number | SNV的数目 |
| INDEL number | INDEL的数目 |
| SNV precision | SNV与标准集比较的precision |
| INDEL precision | INDEL的与标准集比较的precision |
| SNV recall | SNV与标准集比较的recall |
| INDEL recall | INDEL的与标准集比较的recall |

与标准集比较的家庭单元整合结果**precision.recall.txt**

| 列名 | 说明 |
| ----------------- | ------------------------------------------------------------ |
| Family | project.rep.SNV/INDEL；project指的是APP输入的project字段，rep是指第几组家庭 |
| Average Precision | 该组家庭的precision的平均值 |
| Average Recall | 该组家庭的recall的平均值 |
| Precison SD | 该组家庭的precision的SD |
| Recall SD | 该组家庭的recall的SD |

####4 Quartet家系关系评估 ${project}.mendelian.txt

| 列名 | 说明 |
| ----------------------------- | ------------------------------------------------------------ |

+ 17
- 15
codescripts/extract_tables.py Просмотреть файл

@@ -133,21 +133,23 @@ benchmark = benchmark.round(2)
benchmark.to_csv('variants.calling.qc.txt',sep="\t",index=0)

all_rep = [x.split('_')[6] for x in benchmark['Sample']]
rep = list(set(all_rep))
columns = ['Family','Average Precision','Average Recall','Precison SD','Recall SD']
df_ = pd.DataFrame(columns=columns)
for i in rep:
string = "_" + i + "_"
sub_dat = benchmark[benchmark['Sample'].str.contains('_1_')]
mean = list(sub_dat.mean(axis = 0, skipna = True))
sd = list(sub_dat.std(axis = 0, skipna = True))
family_name = project_name + "." + i + ".SNV"
df_ = df_.append({'Family': family_name, 'Average Precision': mean[0], 'Average Recall': mean[2], 'Precison SD': sd[0], 'Recall SD': sd[2] }, ignore_index=True)
family_name = project_name + "." + i + ".INDEL"
df_ = df_.append({'Family': family_name, 'Average Precision': mean[1], 'Average Recall': mean[3], 'Precison SD': sd[1], 'Recall SD': sd[3] }, ignore_index=True)
df_ = df_.round(2)
df_.to_csv('precision.recall.txt',sep="\t",index=0)

if all_rep.count(all_rep[0]) == 4:
rep = list(set(all_rep))
columns = ['Family','Average Precision','Average Recall','Precison SD','Recall SD']
df_ = pd.DataFrame(columns=columns)
for i in rep:
string = "_" + i + "_"
sub_dat = benchmark[benchmark['Sample'].str.contains('_1_')]
mean = list(sub_dat.mean(axis = 0, skipna = True))
sd = list(sub_dat.std(axis = 0, skipna = True))
family_name = project_name + "." + i + ".SNV"
df_ = df_.append({'Family': family_name, 'Average Precision': mean[0], 'Average Recall': mean[2], 'Precison SD': sd[0], 'Recall SD': sd[2] }, ignore_index=True)
family_name = project_name + "." + i + ".INDEL"
df_ = df_.append({'Family': family_name, 'Average Precision': mean[1], 'Average Recall': mean[3], 'Precison SD': sd[1], 'Recall SD': sd[3] }, ignore_index=True)
df_ = df_.round(2)
df_.to_csv('precision.recall.txt',sep="\t",index=0)
else:
pass




Редактирование Предпросмотр
Загрузка…
Отмена
Сохранить