实验记录 | DNA数据的处理流程

目前，尝试对代码文件进行一步一步的拆解。去看他是怎样一步一步的做判断的？

我感觉我的状态有点不对劲，内心的那种动力没有了。反而内心会被一些东西困扰，这些东西阻碍了我的前进。我花半个小时，理一下。
我内心困惑什么呢？还是昨天和老妈讨论的，以及今天和董舜讨论的。
（1）我觉得自己在课题组中被“边缘化”。

没有被纳入到正式群中
做的事情也是细枝末节
目前也没有参与到一些正式的成员所参与的一些事情中，比如组会讲文献，汇报课题
我所看到的和实际上的
所有的同学都渴望加入到这个课题组中，让我产生了一种危机感

这些，让我不太舒服。我害怕重蹈肖老师的覆辙，我很紧张。所以，现在的我，该怎样与这些疙瘩和解呢？
（1）首先，明确一点，你现在不过才过来一个星期。和老师互相都不太熟悉。老师，不会把一些重要的任务交给你是最正常不过的事情。如果是我，我会先安排一些简单的事情，看看她做的怎么样？时间久了之后，如果这个小姑娘我很喜欢，那么我会交给她一些任务。信任的建立是需要时间的，不可能一开始就居于“核心”的位置，要学会忍耐，要学会蓄光。
（2）你现在着手所做的事情虽然是比较基础的，但是是与这个课题的主线相关的一条支线。也非常的重要。将来，也很有可能被纳入到这个体系之中。所以，不存在“游离”一说。你不必有如此的小情绪。
（3）你现在虽然并未是正式的成员，但是是存在可能性的。因为，你的本科有相关的分析的基础，你自己并不差的。老师，主要看两个方面的内容，一方面是你的专业背景与实验室的契合程度，另一方面就是你的专业能力等其它方面的能力。你自己并不差的。
（4）而且，工作的过程中，你的确感觉到快乐不是。不要去在意最终的结果，而且你来到这里的初心就是多学一点东西，你现在的确正在学习中，这难道不是一件快乐的事情吗？值得满足了。
（5）所以，最终，你只需要踏踏实实的安定下来就够了。“请君看取榖纹平”，踏实的过好每一天，每一天都在进步，都比昨天多学习一点东西，就足够了。
以上。

我现在卡在什么地方了呢？
（1）首先，我并不理解为什么四个GATK的结果文件为空？
（2）其次，我并不理解为什么最终没有跑出来mutation calling的结果文件（如例子中展示的）？
这两个问题，并不是纯粹的技术性的问题，而是涉及到整个流程的原理性的问题。需要返回到源代码，去查看它的处理过程（必然经历的一个步骤）。

一个人终究要经过丑陋走向美丽的。

所以，我现在就入手去整理它的每一步的分析过程。
一、准备阶段
（1）定义过程中使用到的参数
（2）准备程序、参考文件
（3）准备工作文件夹

二、比对阶段
非常简洁，只是使用了一个函数alignment（）。

判断输入的文件的类型（这是接下来一切处理步骤的前提）：

exome-seq：rna=0（本次，我们只论述输入文件为dna的处理情况）
rna-seq：rna=1
deep exome-seq：rna=2

如果输入的是bam文件，想将其转换为fastq文件。

序列比对：
/home/xxzhang/workplace/software/bwa/bwa mem -v 1 -t 32 -a -M ./geneome/hg19/hg19.fa ./output/tumor/fastq1.fastq ./output/tumor/fastq2.fastq > ./output/tumor/alignment.sam

mem:使用mem比对算法。
-v:
-a:将所有的比对结果都输出，包括 single-end 和 unpaired paired-end的 reads，但是这些比对的结果会被标记为次优。

（这里过段时间，找到help文档的时候再进行注释）

输入文件：fastq1.fastq；fastq2.fastq
输出文件：alignment.sam

add read group（并不是特别懂）:

/home/xxzhang/workplace/software/java/jdk1.8.0_291/bin/java -jar /home/xxzhang/workplace/QBRC/somatic_script/picard.jar AddOrReplaceReadGroups INPUT=./output/tumor/alignment.sam OUTPUT=./output/tumor/rgAdded.bam SORT_ORDER=coordinate RGID=tumor RGLB=tumor RGPL=illumina RGPU=tumor RGSM=tumor CREATE_INDEX=true VALIDATION_STRINGENCY=LENIENT

输入文件：alignment.sam
输出文件：rgAdded.bam

参考链接：https://broadinstitute.github.io/picard/
这行代码中使用了工具“AddOrReplaceReadGroups”，下面对这些内容进行介绍：
https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups
我不太明白什么是read group?

read group：
Replace read groups in a BAM file.This tool enables the user to replace all read groups in the INPUT file with a single new read group and assign all reads to this read group in the OUTPUT BAM file.

标记重复：
/home/xxzhang/workplace/software/java/jdk1.8.0_291/bin/java -jar /home/xxzhang/workplace/QBRC//somatic_script/picard.jar MarkDuplicates INPUT=./output/tumor/rgAdded.bam OUTPUT=./output/tumor/dupmark.bam CREATE_INDEX=true VALIDATION_STRINGENCY=STRICT REMOVE_SEQUENCING_DUPLICATES=true METRICS_FILE=./output/tumor/dupmark_metrics.txt

输入文件：rgAdded.bam
输出文件：dupmark.bam

这行代码中使用了工具“MarkDuplicates”，这行代码的意思，顾名思义，就是标记重复。我觉得我的错误也不在这个环节。

indel重新比对：

我们的物种是人类，因此使用的比对文件是mills和1000g这两个文件。
比对前：先对使用过的文档排一个顺序。
#/home/xxzhang/workplace/software/java/jdk1.8.0_291/bin/java -Djava.io.tmpdir=./output/tumor/tmp -jar /home/xxzhang/workplace/QBRC//somatic_script/picard.jar ReorderSam INPUT=./output/tumor/dupmark.bam OUTPUT=./output/tumor/dupmark.bam SEQUENCE_DICTIONARY=./geneome/hg19/hg19.dict CREATE_INDEX=TRUE

这行代码其实是可有可无。注释掉，我们真正需要调整顺序的是参考的vcf文件。

RealignerTargetCreator
/home/xxzhang/workplace/software/java/jdk1.8.0_291/bin/java -Djava.io.tmpdir=./output/tumor/tmp -jar /home/xxzhang/workplace/QBRC//somatic_script/GenomeAnalysisTK.jar -T RealignerTargetCreator -R ./geneome/hg19/hg19.fa --num_threads 32 -known ./geneome/hg19/hg19.fa_resource/Mills_and_1000G_gold_standard.indels.hg19.vcf -known ./geneome/hg19/hg19.fa_resource/1000G_phase1.snps.high_confidence.hg19.vcf -o ./output/tumor/tumor_intervals.list -I ./output/tumor/dupmark.bam > ./output/tumor/index.out

其实这一步已经有生成结果文件tumor_intervals.txt，只不过不能够将屏显转移到index.out,并生成它。所以，从整体上看，是没有问题的。

IndelRealigner
/home/xxzhang/workplace/software/java/jdk1.8.0_291/bin/java -Djava.io.tmpdir=./output/tumor/tmp -jar /home/xxzhang/workplace/QBRC//somatic_script/GenomeAnalysisTK.jar -T IndelRealigner --filter_bases_not_stored --disable_auto_index_creation_and_locking_when_reading_rods -R ./geneome/hg19/hg19.fa -known ./geneome/hg19/hg19.fa_resource/Mills_and_1000G_gold_standard.indels.hg19.vcf -known ./geneome/hg19/hg19.fa_resource/1000G_phase1.snps.high_confidence.hg19.vcf -targetIntervals ./output/tumor/tumor_intervals.list -I ./output/tumor/dupmark.bam -o ./output/tumor/realigned.bam >./output/tumor/tumor_realign.out

这一步也生成了结果文件。

base recalibration：

BaseRecalibrator
/home/xxzhang/workplace/software/java/jdk1.8.0_291/bin/java -Djava.io.tmpdir=./output/tumor/tmp -jar /home/xxzhang/workplace/QBRC//somatic_script/GenomeAnalysisTK.jar -T BaseRecalibrator -R ./geneome/hg19/hg19.fa -knownSites ./geneome/hg19/hg19.fa_resource/dbsnp.hg19.vcf -knownSites ./geneome/hg19/hg19.fa_resource/Mills_and_1000G_gold_standard.indels.hg19.vcf -I ./output/tumor/realigned.bam -o ./output/tumor/tumor_bqsr > ./output/tumor/table.out
PrintReads
/home/xxzhang/workplace/software/java/jdk1.8.0_291/bin/java -Djava.io.tmpdir=./output/tumor/tmp -jar /home/xxzhang/workplace/QBRC//somatic_script/GenomeAnalysisTK.jar -T PrintReads -rf NotPrimaryAlignment -R ./geneome/hg19/hg19.fa -I ./output/tumor/realigned.bam -BQSR ./output/tumor/tumor_bqsr -o ./output/tumor/tumor.bam > ./output/tumor/tumor_recal.out

这两步都没有问题。

进一步处理bam文件：
/home/xxzhang/workplace/software/sambamba/sambamba index -t 32 ./output/tumor/tumor.bam
/home/xxzhang/miniconda3/bin/bcftools mpileup -I -f ./geneome/hg19/hg19.fa ./output/tumor/tumor.bam > ./output/tumor/tumor.mpileup

三、突变筛选阶段
/home/xxzhang/workplace/software/strelka/bin/configureStrelkaGermlineWorkflow.py --bam /home/xxzhang/workplace/QBRC/output/tumor/tumor.bam --referenceFasta ./geneome/hg19/hg19.fa --runDir ./output/strelka/ --exome

输入文件：tumor.bam（也是上一步生成的）
输出文件：

/home/xxzhang/workplace/QBRC/output/strelka/runWorkflow.py -m local -j 32

这一步觉得也有问题，为什么只是calling出了一行？对，我觉得问题也在这里。

四、整合vcf文件

过滤vcf文件

/home/xxzhang/workplace/software/R/R-3.6.3/bin/Rscript /home/xxzhang/workplace/QBRC//somatic_script/filter_vcf.R ./output NA

annovar注释

/home/xxzhang/workplace/QBRC/annovar/table_annovar.pl ./output/germline_mutations.txt /home/xxzhang/workplace/QBRC/annovar//humandb/ -buildver hg19 -out ./output/germline_mutations -remove -protocol refGene,ljb26_all,cosmic70,esp6500siv2_all,exac03,1000g2015aug_all -operation g,f,f,f,f,f -nastring

/home/xxzhang/workplace/software/R/R-3.6.3/bin/Rscript /home/xxzhang/workplace/QBRC//somatic_script/filter_vcf2.R ./output/somatic_mutations.txt ./output/somatic_mutations.hg19_multianno.txt ./output/somatic_mutations_hg19.txt somatic

/home/xxzhang/workplace/QBRC/annovar/annotate_variation.pl -geneanno -dbtype refGene -buildver hg19 ./output/somatic_mutations_hg19.txt /home/xxzhang/workplace/QBRC/annovar//humandb/

/home/xxzhang/workplace/QBRC/annovar/coding_change.pl --includesnp --alltranscript --newevf ./output/somatic_mutations_hg19.txt_tmp.txt ./output/somatic_mutations_hg19.txt.exonic_variant_function /home/xxzhang/workplace/QBRC/annovar//humandb//hg19_refGene.txt /home/xxzhang/workplace/QBRC/annovar//humandb//hg19_refGeneMrna.fa >/dev/null 2>/dev/null

/home/xxzhang/workplace/software/R/R-3.6.3/bin/Rscript /home/xxzhang/workplace/QBRC//somatic_script/add_fs_annotation.R ./output hg19 somatic

注释内容：用到了什么文件？

实验记录 | DNA数据的处理流程相关推荐

实验记录 | scATAC-seq数据的比对（一）
1.熟悉scATAC-seq的数据的基本格式 ATAC-seq数据,用于测序的reads如图所示. 按照它的参考文件中的描述,一共有4条测序的reads. (1)双端测序的Read 1N和Read 2 ...
《保护我们的数字遗产：DNA数据存储》白皮书发布
编者按: 2020年10月,Twist Bioscience.Illumina.Western Digital(西部数据).微软研究院等公司和机构联合成立DNA数据存储联盟(DNA data stor ...
HIT-大数据分析Lab1：数据预处理-实验记录
本文是哈工大大数据分析实验1的完整实验记录,包含环境配置,相关知识介绍以及实验解析.希望对后来人有帮助(新手小白没什么头绪,走一步查一步对应的博客o(╥﹏╥)o),博客链接之间会穿插一些我自己的理解, ...
【实验记录】EA-MLP(演化算法--全连接神经网络)实验记录
large scale evaluation net -- MLP全连接实验记录 Ⅰ. Experiment detail Ⅱ. Method Vertex Edge DNA Evolution_po ...
CSAPP Lab2 实验记录 ---- Bomb Lab（Phase 1 - Phase 6详细解答 + Secret Phase彩蛋解析）
文章目录 Lab 总结博客链接实验前提引子实验需要指令及准备 Phase 1 Phase 2 Phase 3 Phase 4 Phase 5 Phase 6 Phase Secret(彩蛋Phas ...
操作系统真象还原实验记录之实验十一：实现中断处理（二）
操作系统真象还原实验记录之实验十一:实现中断处理(二) 书p335 7.6.2 改进中断处理程序,并调快时钟 1.实验代码第一次修改对应书p335 7.6.2 改进中断处理程序这次是上一次实验的 ...
用Python实现不同数据源的对象匹配【实验记录】
任务简介: 现有两份针对同一主题的数据,但是在人物的属性名称及格式上有所不同,需要对两份数据进行匹配来确定是同一个人. 匹配属性: 1 人名 2 出生日期 3 国籍原始数据举例 1 数据源1(以下简 ...
中国电子实验记录（ELN）系统行业研究报告（2022版）
内容简介: 电子实验记录(ELN)系统是伴随着无纸化管理在实验室的应用需求而逐步发展起来的专业化系统开发平台.该系统以人为中心,可以专门用来进行辅助实验设计并结构化地记录实验过程与数据,支持实验室研发 ...
论文推荐：SoLid 通过让公民控制自己的数据简化政府流程
SoLiD 是一个令人兴奋的新项目,由万维网发明者 Tim Berners-Lee 爵士在麻省理工学院启动. 该项目旨在从根本上改变 Web 应用程序的中心化趋势, 它将真正地让数据所有权属于用户,并 ...

实验记录 | DNA数据的处理流程

实验记录 | DNA数据的处理流程相关推荐

最新文章

热门文章