实验记录 |somatic.pl的运行2

一个同学分享了一个网盘文件，上面有gatk相关的文件，所以我尝试下载，是否能够使用。
参考链接：https://blog.csdn.net/xxxie_/article/details/100111991
我的心态是崩的，原来早就有同学分享到网盘里了。
我也有匆匆浏览过这个网页，但是没有注意到它在文后的一个分享。
还是要静下来去做事吧！

现在已经有一个文件了，把它放在合适的文件夹下，重新跑程序，看看结果。
移动完成。
重新跑代码。
perl somatic.pl NA NA example/example_dataset/sequencing/SRR7246238_1.fastq.gz example/example_dataset/sequencing/SRR7246238_2.fastq.gz 32 hg19 /home/zxx/QBRC/geneome/hg19/hg19.fa /usr/lib/jvm/jdk1.8.0_181/bin/java human 1 /home/zxx/QBRC/disambiguate_pipeline

毫无疑问，又出现了一堆的错误，但是希望没有之前的问题。

MESSAGE: Could not read file /home/zxx/QBRC/geneome/hg19/hg19.fa_resource/1000G_phase1.snps.high_confidence.hg19.vcf because file ‘/home/zxx/QBRC/geneome/hg19/hg19.fa_resource/1000G_phase1.snps.high_confidence.hg19.vcf’ does not exist

ERROR MESSAGE: Could not read file /home/zxx/QBRC/geneome/hg19/hg19.fa_resource/dbsnp.hg19.vcf because file ‘/home/zxx/QBRC/geneome/hg19/hg19.fa_resource/dbsnp.hg19.vcf’ does not exist

ERROR MESSAGE: Could not read file /home/zxx/QBRC/human/tumor/tumor_bqsr because it does not exist

还是配置文件的问题，难道不能一起出现吗？
好的，我继续下载这一部分。

下载完成（2021/05/15）。
现在将移动硬盘中的文件移动到指定目录下。
经过上次的经验积累，我发现用正规途径（移动指令或交互式界面操作）比手动的拖放，更不容易造成文件的损坏。
所以，我现在用文件的移动指令，在进行此次操作。
类如：
cp dbsnp_138.hg19.vcf.idx.gz /home/zxx/QBRC/geneome/hg19/hg19.fa_resource
这样操作之后的话，就是在指定的文件夹下进行解压。我学会了批量的操作。
因为是.gz结尾的文件，于是选择使用解压指令gunzip。
gunzip *.gz
我很喜欢这种批量的操作。
另外，我最近还想学习一下shell语言中的循环是怎么一回事。感觉直接用shell脚步批量跑很方便。之前有接触过，但是没有认真的研究过。
截止到现在，我的四个文件复制移动完成。
我决定重新运行一下代码。
不过，在此之前，需要改一下dbsnp文档的名字。
mv dbsnp_138.hg19.vcf dbsnp.hg19.vcf
后来报错显示，1000那个文档也需要改名字——粗心的我。
继续改名字。
mv 1000G_phase1.snps.high_confidence.hg19.sites.vcf 1000G_phase1.snps.high_confidence.hg19.vcf
实践证明，还是交互式的处理比较方便（更加适合现在的我）。

好的，现在切换到somatic.pl所在的路径下，重新运行代码。
perl somatic.pl NA NA example/example_dataset/sequencing/SRR7246238_1.fastq.gz example/example_dataset/sequencing/SRR7246238_2.fastq.gz 32 hg19 /home/zxx/QBRC/geneome/hg19/hg19.fa /usr/lib/jvm/jdk1.8.0_181/bin/java human 1 /home/zxx/QBRC/disambiguate_pipeline

另外，参考昨天那个somatic.pl的calling流程，如果是gatk相关文件缺失导致报错的话，说明程序已经运行到中间步骤了。
重新运行之后，产生了新的错误。但幸运的是，不是以前的错误了。
我重新进行分析。

ERROR MESSAGE: Lexicographically sorted human genome sequence detected in reference. Please see http://gatkforums.broadinstitute.org/discussion/58/companion-utilities-reordersamfor more information.

Error details: reference contigs = [chr10, chr11, chr11_gl000202_random, chr12, chr13, chr14, chr15, chr16, chr17_ctg5_hap1, chr17, chr17_gl000203_random, chr17_gl000204_random, chr17_gl000205_random, chr17_gl000206_random, chr18, chr18_gl000207_random, chr19, chr19_gl000208_random, chr19_gl000209_random, chr1, chr1_gl000191_random, chr1_gl000192_random, chr20, chr21, chr21_gl000210_random, chr22, chr2, chr3, chr4_ctg9_hap1, chr4, chr4_gl000193_random, chr4_gl000194_random, chr5, chr6_apd_hap1, chr6_cox_hap2, chr6_dbb_hap3, chr6, chr6_mann_hap4, chr6_mcf_hap5, chr6_qbl_hap6, chr6_ssto_hap7, chr7, chr7_gl000195_random, chr8, chr8_gl000196_random, chr8_gl000197_random, chr9, chr9_gl000198_random, chr9_gl000199_random, chr9_gl000200_random, chr9_gl000201_random, chrM, chrUn_gl000211, chrUn_gl000212, chrUn_gl000213, chrUn_gl000214, chrUn_gl000215, chrUn_gl000216, chrUn_gl000217, chrUn_gl000218, chrUn_gl000219, chrUn_gl000220, chrUn_gl000221, chrUn_gl000222, chrUn_gl000223, chrUn_gl000224, chrUn_gl000225, chrUn_gl000226, chrUn_gl000227, chrUn_gl000228, chrUn_gl000229, chrUn_gl000230, chrUn_gl000231, chrUn_gl000232, chrUn_gl000233, chrUn_gl000234, chrUn_gl000235, chrUn_gl000236, chrUn_gl000237, chrUn_gl000238, chrUn_gl000239, chrUn_gl000240, chrUn_gl000241, chrUn_gl000242, chrUn_gl000243, chrUn_gl000244, chrUn_gl000245, chrUn_gl000246, chrUn_gl000247, chrUn_gl000248, chrUn_gl000249, chrX, chrY

ERROR MESSAGE: Could not read file /home/zxx/QBRC/human/tumor/realigned.bam because java.io.FileNotFoundException: human/tumor/realigned.bam (No such file or directory)

ERROR MESSAGE: Could not read file /home/zxx/QBRC/human/tumor/tumor_bqsr because it does not exist

这个程序是串联的，所以，后面的错误往往是前面的错误所引起的，我现在主要要做的就是解决前面出现的这个错误。

具体的解决思路：
百度翻译了一下，报错的主要的意思是：

参考基因组中检测到按字典排序的人类基因组序列

了解到这一点之后，我们看一下具体的细节。果然是字典排序的方式。

chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr1, chr20, chr21, chr22, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chrM, chrUn,chrX, chrY

这里补充一下，什么是字典排序（这个之前建索引的时候有接触过）？
参考链接：百度百科-字典排序

设想一本英语字典里的单词，何者在前何者在后？
显然的做法是先按照第一个字母、以 a、b、c……z 的顺序排列；如果第一个字母一样，那么比较第二个、第三个乃至后面的字母。如果比到最后两个单词不一样长（比如，sigh 和 sight），那么把短者排在前。
通过这种方法，我们可以给本来不相关的单词强行规定出一个顺序。“单词”可以看作是“字母”的字符串，而把这一点推而广之就可以认为是给对应位置元素所属集合分别相同的各个有序多元组规定顺序。

所以，在这里的话，是不太希望用字典排序这种方式吗？
我预估问题出现在了索引上。
我们打开看一下我们上次建的两个参考基因组的索引（dict与fai）。
head hg19.dict

@HD VN:1.6
@SQ SN:chr10 LN:135534747 M5:988c28e000e84c26d552359af1ea2e1d UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr11 LN:135006516 M5:98c59049a2df285c76ffb1c6db8f8b96 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr11_gl000202_random LN:40103 M5:06cbf126247d89664a4faebad130fe9c UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr12 LN:133851895 M5:51851ac0e1a115847ad36449b0015864 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr13 LN:115169878 M5:283f8d7892baa81b510a015719ca7b0b UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr14 LN:107349540 M5:98f3cae32b2a2e9524bc19813927542e UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr15 LN:102531392 M5:e5645a794a8238215b2cd77acb95a078 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr16 LN:90354753 M5:fc9b1a7b42b97a864f56b348b06095e6 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr17_ctg5_hap1 LN:1680828 M5:d89517b400226d3b56e753972a7ca

head hg19.fa.fai

chr10 135534747 7 50 51
chr11 135006516 138245456 50 51
chr11_gl000202_random 40103 275952126 50 51
chr12 133851895 275993039 50 51
chr13 115169878 412521979 50 51
chr14 107349540 529995262 50 51
chr15 102531392 639491800 50 51
chr16 90354753 744073827 50 51
chr17_ctg5_hap1 1680828 836235693 50 51
chr17 81195210 837950145 50 51

从上面两个文件的内容看，确实都是报错中显示的“字典排序”的类型。但是，除了字典排序，还有什么排序类型呢？
而且我建立参考索引的过程是按照官网进行操作的。
那么，到这里，我们还有最后一点线索，就是错误中提示的去官网： http://gatkforums.broadinstitute.org/discussion/ ，查看是否有比较好的建议。
以“Lexicographically”作为关键词进行搜索，找到了5篇与之相关的文章。

其中，可以参考一下这篇文章。或许可以帮助我们更好的理解运行的机制（2021/05/15 9：20）。
https://gatk.broadinstitute.org/hc/en-us/articles/360057440671-SortSam-Picard-

一般情况下，参考基因组的排序方式有两种：
（1）第1种是按照mapping到的参考基因组的坐标上下游顺序来排序，是samtools的默认排序方法
（2）第2种是按照reads name进行排序，需要增加一个-n参数。

我们来分别尝试体验一下。上一次我们生成索引，采用的指令是：
samtools faidx hg19.fa
这样一种构建索引的方式，是如何构建的呢？
这个索引好像是直接通过samtools对fasta进行排序的。而fasta又来源于我们原始的那些染色体序列，通过cat合并而成。

所以，究竟怎样才好呢？我在这个地方卡住了。

48265987+Edison19991109@users.noreply.github.com
33228841+tianshilu@users.noreply.github.com
curl “https://api.github.com/users/MenglinC/repos?type=owner&sort=updated” -s
| sed -nE ‘s#^.“name”: “([^”]+)",.KaTeX parse error: Expected 'EOF', got '#' at position 1: #̲\1#p' \ …#\2#p’|sort -u
如何通过github用户名找到该用户的邮箱
https://api.github.com/users/Edison19991109/events/public

经过上面一段尝试，我决定使用VPN。
于是联系了某宝客服，在12:08，我终于自己能够人生中第一次连上谷歌了。
我激动的心情，就像是村里通电了一样。
不得不说，好的工具能够让我们便宜许多。
以“Lexicographically sorted human genome sequence detected in reference.”为关键词进行搜索：
得到如下解答：
可以同时比较一下bing的对于同一个问题的检索结果。所以总结得出一个好的搜索引擎是多么的重要。
所以，我买的是一周的时间。在这一周的时间内，要好好利用资源。
但是，比较遗憾的是这个vpn只能应用于window平台。
所以，我只能在window这边找到比较合适的解决方案之后，在Linux那边集中的处理与尝试。

这个问题，具体的建议如下：
参考链接：https://www.biostars.org/p/63337/

1. MESSAGE: Lexicographically sorted human genome sequence detected in reads.
   For safety’s sake the GATK requires human contigs in karyotypic order: 1, 2, …, 10, 11, …, 20, 21, 22, X, Y with M either leading or trailing these contigs

所以，正常的顺序应该是：

chr1,chr2,chr3,chr4 etc.

参考链接：https://www.biostars.org/p/331664/
这篇文章就在探讨，如何将自己的数据按照这种方式进行排序。
参考链接：https://gatk.broadinstitute.org/hc/en-us/articles/360035532232?id=1328

Version-specific alert for GATK 3.5
In version 3.5, we added some beefed-up VCF sequence dictionary validation. Unfortunately, as a side effect of the additional checks, some users have experienced an error that starts with “ERROR MESSAGE: Lexicographically sorted human genome sequence detected in variant.” that is due to unintentional activation of a check that is not necessary. This will be fixed in the next release; in the meantime -U ALLOW_SEQ_DICT_INCOMPATIBILITY can be used (with caution) to override the check.

重新sort我们的fasta格式的文件（我觉得这就是我要找的答案）。

• Sort your reference

samtools faidx reference.fa

`cut -f1 reference.fa.fai|sort -k1V|parallel -k 'samtools faidx reference.fa {}' > ref_sorted.fa`

• Create a dictionary for the sorted ref
  java -jar picard.jar CreateSequenceDictionary R=ref_sorted.fa O=ref_sorted.dict

• Reorder bam file
  java -jar picard.jar ReorderSam I=original.bam O=reordered.bam R=ref_sorted.fa CREATE_INDEX=TRUE

选择尝试这种解决方案（记得备份），觉得归根到底还是自己的基础太差，不会使用指令进行文件操作，等一下我要好好分析一下上面这个shell指令的意思。

首先cut指令在shell中，是截取字符串的意思。参数-f1截取第一个字段。对应到文件中，就是字符chr10等等。

chr10 135534747 7 50 51
chr11 135006516 138245456 50 51
chr11_gl000202_random 40103 275952126 50 51
chr12 133851895 275993039 50 51
chr13 115169878 412521979 50 51

|在linux中的意思是分割两个命令的管道符。
sort 是linux中的排序指令。-k1指定为第一列。
parallel多进程并行运算。
修改代码为：
cut -f1 hg19.fa.fai | sort -k1V|parallel -k 'samtools faidx hg19.fa {}' > hg19_sorted.fa

得到了文件hg19_sorted.fa并将其更名为hg19.fa。
mv hg19_sorted.fa hg19.fa

在这个基础上，重新创建索引文件。

• hg19.dict
gatk CreateSequenceDictionary -R hg19.fa
  需要一定的时间。创建完成之后，我们看一下，结果文件是否是我们想要的顺序。
  有点开心，是我们想要的这种格式。

@HD VN:1.6
@SQ SN:chr1 LN:249250621 M5:1b22b98cdeb4a9304cb5d48026a85128 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr1_gl000191_random LN:106433 M5:d75b436f50a8214ee9c2a51d30b2c2cc UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr1_gl000192_random LN:547496 M5:325ba9e808f669dfeee210fdd7b470ac UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr2 LN:243199373 M5:a0d9851da00400dec1098a9255ac712e UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr3 LN:198022430 M5:641e4338fa8d52a5b781bd2a2c08d3c3 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr4 LN:191154276 M5:23dccd106897542ad87d2765d28a19a1 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr4_ctg9_hap1 LN:590426 M5:fa24f81b680df26bcfb6d69b784fbe36 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr4_gl000193_random LN:189789 M5:dbb6e8ece0b5de29da56601613007c2a UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr4_gl000194_random LN:191469 M5:6ac8f815bf8e845bb3031b73f812c012 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr5 LN:180915260 M5:0740173db9ffd264d728f32784845cd7 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr6 LN:171115067 M5:1d3a93a248d92a729ee764823acbbc6b UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa
@SQ SN:chr6_apd_hap1 LN:4622290 M5:fe71bc63420d666884f37a3ad79f3317 UR:file:/home/zxx/QBRC/geneome/hg19/hg19.fa

以同样的方式，依据新的参考基因组，创建fai参考索引。

• hg19.fa.fai
samtools faidx hg19.fa
  同样的，也是我们想要的那种格式。

chr1 249250621 6 60 61
chr1_gl000191_random 106433 253404827 60 61
chr1_gl000192_random 547496 253513056 60 61
chr2 243199373 254069683 60 61
chr3 198022430 501322385 60 61
chr4 191154276 702645195 60 61
chr4_ctg9_hap1 590426 896985392 60 61
chr4_gl000193_random 189789 897585681 60 61
chr4_gl000194_random 191469 897778656 60 61
chr5 180915260 897973323 60 61
chr6 171115067 1081903844 60 61
chr6_apd_hap1 4622290 1255870844 60 61
chr6_cox_hap2 4795371 1260570188 60 61
chr6_dbb_hap3 4610396 1265445497 60 61

不得不说GATK相当挑剔。

所以，到这里，再次运行命令行。
perl somatic.pl NA NA example/example_dataset/sequencing/SRR7246238_1.fastq.gz example/example_dataset/sequencing/SRR7246238_2.fastq.gz 32 hg19 /home/zxx/QBRC/geneome/hg19/hg19.fa /usr/lib/jvm/jdk1.8.0_181/bin/java human 1 /home/zxx/QBRC/disambiguate_pipeline

当然，又出了一堆的问题。
剩下的报错，放到下一个实验记录里解释。