实验记录 | 8/3

今天早上目标：

（1）将OX1931、CML656两个病人的sra文件转换为fastq文件。
（2）根据注释信息，选择我们想要进一步测试的来源于不同患者的fastq文件。

for i in $(cat ls.log); do  fastq-dump --split-files --gzip  /home/xxzhang/workplace/QBRC/data/CML/OX1931/$i/$i.sra; done 
for i in $(cat ls.log); do  fastq-dump --split-files --gzip  /home/xxzhang/workplace/QBRC/data/CML/CML656/$i/$i.sra; done 

系统显示说，没有ls.log的文件。

cat: ls.log: No such file or directory

让我来仔细看一下，我之前用的ls.log文件的内容是什么？

SRR3052548
SRR3052549
SRR3052550
SRR3052551
SRR3052552
SRR3052553
……

因为好久没有处理这件事了。所以，有一些文件的内容基本上都忘记了。
本来想利用以下指令，重新生成。但是发现我的目录下，已经有存储这个信息的文件了。

ls >ls.log 

想想，居然都是一个月前的事情了。时间过得真快，时间也确实可以改变一些东西，希望自己能够利用这些时间去做有意义的事情，去解决有意思的问题。老师说这个流程不好做，那我偏要弄下去。我相信是可以有结果的，实在不行，拿来练手也可以。反正我是不会放弃的。

出现了问题：

2021-08-03T02:43:16 fastq-dump.2.11.0 err: error unexpected while resolving query within virtual file system module - No accession to process ( 500 )
Failed to call external services.
2021-08-03T02:43:17 fastq-dump.2.11.0 err: error unexpected while resolving query within virtual file system module - No accession to process ( 500 )
Failed to call external services.

问题出在了这里。

for i in $(cat CML656_cellset.txt); do  fastq-dump --split-files --gzip  /home/xxzhang/workplace/QBRC/data/CML/CML656/$i/$i.sra; done

这个批量的处理，并没有成功。我想了一下，也许我其实并不需要全部把它解压下来，我只需要从中挑几个来做就可以了。

---

所以，现在的任务变成了，从中挑几个来处理。
由于两个数据文件的信息并不是特别全，所以现在可能需要写个code，将他们一一对应起来。merge的标准是GSM号的对应关系。

在注释文件中，选择了我比较感兴趣的具有生物学差异的样本，整理成了下图的list。可以有四个维度的比较（作为一个强迫症，非要证实一下不可，看看从突变这个维度看，能不能找出点啥有趣的东西）。

---

确定这个list之后，使用sratoolkit进行下载。
将这个sra文件，整理成一个list，然后批量下载。

prefetch --option-file /home/xxzhang/workplace/QBRC/data/CML/TESTclass/SRR_Acc_List.txt --output-directory /home/xxzhang/workplace/QBRC/data/CML/TESTclass

正在下载中。
下载完成（12:16）。

还是想试一下，批量的fastq-dump行不行？

for i in $(cat ls.log); 
do  fastq-dump --split-files --gzip  /home/xxzhang/workplace/QBRC/data/CML/TESTclass/$i/$i.sra; 
done 

成功了，我觉得之前错误的主要原因是由于我没有将它的格式整理一下。

(base) [xxzhang@mu02 TESTclass]$ 
for i in $(cat ls.log);
> do  fastq-dump --split-files --gzip  /home/xxzhang/workplace/QBRC/data/CML/TESTclass/$i/$i.sra;
> done

（14：00）早上的目标到这里基本上算是完成了，基本上比较顺利。

---

现在制定下午和晚上的目标。下午（16：00）要开组会，要把这部分的时间空出来。
（1）首先核实perl代码的每一个参数（因为我之前测试参数的时候有修改过，但是现在隔了那么久都忘记了，为了防止出错，我先用一个进行测试）
（2）另外一方面，我想对这个代码文件进行修改（我现在比较想做的一个工作，是对这个pipeline进行优化，虽然现在任务难度降低，但是我们可以自己在这方面做些有意思的事情）

• 每一个环节标注出来时间，到时候，将这部分的运行时间的数据和我们的数据的大小进行比较（进行一些比较有意思的探索）。
• 保留比较重要的过程性的文件（这个pipeline可能介于存储空间的考虑，所以把过程性文件删掉了，但是咱们想要知道，是否有必要）。但是这个具体要怎么做？我需要多看看文献。

Li H. Toward better understanding of artifacts in variant calling from highcoverage samples. Bioinformatics. 2014;30(20):2843–51. #找到了李恒大佬在n多年写的一篇文章。比较有意思，可以多学习学习别人的优秀的思维。#这些都是免费的优秀的资源，我觉得我现在之所以迷茫是好的东西见的太少。

好了，小菜鸟继续干活。

好了，现在和作者更新过后的流程人为的对照了一下，又对somatic.pl文件做了调整。其中，最大的变动就是：

 #这个我还真没有细看，在他更新后的文件中下面的几行好像是删去了
  #我确认一下是否删去，然后再看一下这几行的用途是什么？#以及它的output是什么?
  
  #system_call("/home/xxzhang/workplace/QBRC/annovar/annotate_variation.pl -geneanno -dbtype refGene -buildver ".$build." ".$output."/".$type."_mutations_".$build.".txt ".$annovar_path.$annovar_db);
  #system_call("/home/xxzhang/workplace/QBRC/annovar/coding_change.pl --includesnp --alltranscript --newevf ".$output."/".$type."_mutations_".$build.".txt_tmp.txt ".$output."/".$type."_mutations_".$build.".txt".
 #   ".exonic_variant_function ".$annovar_path.$annovar_db."/".$build."_refGene.txt ".$annovar_path.$annovar_db."/".$build."_refGeneMrna.fa >/dev/null 2>/dev/null");
 # system_call("/home/xxzhang/workplace/software/R/R-3.6.3/bin/Rscript ".$path."/somatic_script/add_fs_annotation.R ".$output." ".$build." ".$type);
  #system_call("rm -f ".$output."/".$type."_mutations_".$build.".txt?*");

这个先暂时的放一下，其实要弄明白也不难。

somatic.pl文件准备好之后，然后写perl批处理的代码，到时候直接bash即可。

anno<-read.table("cml_table.txt",header = T)
dim(anno)[1]
line1<-rep("perl /home/xxzhang/workplace/QBRC//somatic.pl NA	NA RNA:/home/xxzhang/workplace/QBRC/data/CML/TESTclass/",28)
line2<-rep("_1.fastq.gz	NA 32	hg38	/home/xxzhang/workplace/QBRC/geneome/hg38/hg38.fasta /home/xxzhang/miniconda3/bin/java /home/xxzhang/workplace/QBRC/output/",28)
line3<-rep("/	human	1	./disambiguate_pipeline",28)
data<-cbind(line1,anno$Run,line2,anno$Title,line3)
write.table(data,"test.sh",sep = "",quote = F,col.names = F,row.names = F)

最后得到test.sh的指令。

perl /home/xxzhang/workplace/QBRC//somatic.pl NA	NA RNA:/home/xxzhang/workplace/QBRC/data/CML/TESTclass/SRR3052202_1.fastq.gz	NA 32	hg38	/home/xxzhang/workplace/QBRC/geneome/hg38/hg38.fasta /home/xxzhang/miniconda3/bin/java /home/xxzhang/workplace/QBRC/output/NBMA3/	human	1	./disambiguate_pipeline
perl /home/xxzhang/workplace/QBRC//somatic.pl NA	NA RNA:/home/xxzhang/workplace/QBRC/data/CML/TESTclass/SRR3052203_1.fastq.gz	NA 32	hg38	/home/xxzhang/workplace/QBRC/geneome/hg38/hg38.fasta /home/xxzhang/miniconda3/bin/java /home/xxzhang/workplace/QBRC/output/NBMA5/	human	1	./disambiguate_pipeline
perl /home/xxzhang/workplace/QBRC//somatic.pl NA	NA RNA:/home/xxzhang/workplace/QBRC/data/CML/TESTclass/SRR3052204_1.fastq.gz	NA 32	hg38	/home/xxzhang/workplace/QBRC/geneome/hg38/hg38.fasta /home/xxzhang/miniconda3/bin/java /home/xxzhang/workplace/QBRC/output/NBMB5/	human	1	./disambiguate_pipeline
perl /home/xxzhang/workplace/QBRC//somatic.pl NA	NA RNA:/home/xxzhang/workplace/QBRC/data/CML/TESTclass/SRR3052205_1.fastq.gz	NA 32	hg38	/home/xxzhang/workplace/QBRC/geneome/hg38/hg38.fasta /home/xxzhang/miniconda3/bin/java /home/xxzhang/workplace/QBRC/output/NBMC5/	human	1	./disambiguate_pipeline
……

最后，运行下方的指令即可。

bash tesh.sh