实验记录 | 6/29

使用scSplitter拆分10X 的数据的时候，出现了问题。

EXITING because of FATAL ERROR in reads input: quality string length is not equal to sequence length
@ST-K00126:608:HWNLJBBXX:6:2115:23480:3459
@ST-K00126:608:HWNLJBBXX:6:1212:4005:42513 2:N:0:GATCTCAG
SOLUTION: fix your fastq file
@ST-K00126:608:HWNLJBBXX:7:2115:21592:44816
@ST-K00126:608:HWNLJBBXX:7:1215:3478:11460 2:N:0:GATCTCAG

这个并不是这个Python程序自己报的错，而是STAR的错误。猜想主要的原因可能是：fastq文件出了错。
我们先看一下其比对所使用的的fastq文件。

STAR --runThreadN 20 --outSAMmultNmax 1 --outSAMunmapped Within --outSAMorder PairedKeepInputOrder --genomeDir /home/xxzhang/workplace/QBRC/geneome/hg38/STAR --readFilesIn 1_R5_chunk0000,1_R5_chunk0001,1_R5_chunk0002,1_R5_chunk0003 1_R2_chunk0000,1_R2_chunk0001,1_R2_chunk0002,1_R2_chunk0003     --outFileNamePrefix
/home/xxzhang/workplace/QBRC/data/10X/pbmc2_v1_cells/sams/lane_1

所以，通过上述的代码，我们可以看到，我们要处理的文件是，1_R5_chunk0000这些。
我们现在的问题，并不是说这个问题是否存在，而是如何解决这个问题。相信一句话，欲速则不达。还是应该要将心绪平稳下来。

这个错误的意思就是说，尝试着修复我们的文件。但是，如何修复呢？这就让人觉得很奇怪。
我继续看一下，源代码是如何将fastq拆分的？

---

陷入了瓶颈。现在开始着手处理，关于job_somatic.pl的批处理的代码。
作者在文章中是这样解释的：

perl job_somatic.pl design.txt example_file thread build index java17 disambiguate_pipeline

关于，design.txt需要包括这几行参数：
(1) 正常样本（双端1 或 NA）正常样本2（双端2 或 NA） 肿瘤样本（双端1）肿瘤样本（双端2） 输出文件夹 human

~/seq/1799-01N.R1.fastq.gz ~/seq/1799-01N.R2.fastq.gz ~/seq/1799-01T.R1.fastq.gz ~/seq/1799-01T.R2.fastq.gz ~/out/1799-01/ human
~/seq/1799-02N.R1.fastq.gz ~/seq/1799-02N.R2.fastq.gz ~/seq/1799-02T.R1.fastq.gz ~/seq/1799-02T.R2.fastq.gz ~/out/1799-02/ human
~/seq/1799-03N.R1.fastq.gz ~/seq/1799-03N.R2.fastq.gz ~/seq/1799-03T.R1.fastq.gz ~/seq/1799-03T.R2.fastq.gz ~/out/1799-03/ human

这个文件根据我们自己的实际情况修改即可。

还有一个文件就是example_file。这个文件当初我也是向作者要到的。
具体的内容如下：

#!/bin/bash
#ATCH --job-name=c
#BATCH --partition=256GB
#SBATCH --nodes=1
#SBATCH --time=60-00:00:00
#SBATCH --output=./sbatch_output_%j
#SBATCH --error=./sbatch_error_%j
source ~/.bash_profile
##########JOBSTART###########################
perl job_somatic.pl /project/job_somatic.txt example.sh 32 hg38 /project/data/hg38/hs38d1.fa /cm/shared/apps/java/oracle/jdk1.7.0_51/bin/java 3

关于，这个example.sh的内容也是非常的有意思，虽然并不是特别的清楚文件之间的相互的嵌套关系。
其实我觉得关于这个example.sh我觉得可以只有一行指令，

##########JOBSTART###########################

剩下的内容，都是程序之后写入的。之所以会这样说，是因为在job_somatic.pl的代码文件中，是这样说的：

# write header
  open(HEADER,$example) or die "Cannot find the example shell script!\n";
  while ($line1=<HEADER>)
  {
    if ($line1=~/JOBSTART/) {last;}
    print SCRIPT $line1;
  }
  close(HEADER);
}

以上是我的推断。接下来，我们可以使用mtDNA的那个文件进行测试。
我可以具体再看一下，这个文件中，作者是如何提交参数的？或者说提交的参数，分别都进行了怎样的处理。
总体而言的示例代码：

perl /home2/twang6/software/cancer/somatic/job_somatic.pl \
design.txt \
/project/bioinformatics/Xiao_lab/shared/neoantigen/code/somatic/example/example.sh \
32 hg38 \
/project/shared/xiao_wang/data/hg38/hs38d1.fa \
cm/shared/apps/java/oracle/jdk1.7.0_51/bin/java \
0 2
/project/shared/xiao_wang/software/disambiguate_pipeline

修改为我的示例代码：

perl job_somatic.pl ./somatic_script/design.txt ./somatic_script/example.sh  32 hg19 /home/xxzhang/workplace/QBRC/geneome/hg19/hg19.fa /home/xxzhang/workplace/software/java/jdk1.7.0_80/bin/java 1 3 /home/xxzhang/workplace/QBRC/disambiguate_pipeline

design.txt的主要内容（注意用tab键分隔）：

NA NA /home/xxzhang/workplace/QBRC/example/example_dataset/sequencing/SRR7246238_1.fastq.gz  /home/xxzhang/workplace/QBRC/example/example_dataset/sequencing/SRR7246238_2.fastq.gz  ./output_job/	human
NA NA /home/xxzhang/workplace/QBRC/example/example_dataset/sequencing/SRR7246238_1.fastq.gz  /home/xxzhang/workplace/QBRC/example/example_dataset/sequencing/SRR7246238_2.fastq.gz  ./output_job/	human
NA NA /home/xxzhang/workplace/QBRC/example/example_dataset/sequencing/SRR7246238_1.fastq.gz  /home/xxzhang/workplace/QBRC/example/example_dataset/sequencing/SRR7246238_2.fastq.gz  ./output_job/	human

我觉得example.sh这个文件的可能结果是：

#!/bin/bash
#ATCH --job-name=c
#BATCH --partition=256GB
#SBATCH --nodes=1
#SBATCH --time=60-00:00:00
#SBATCH --output=./sbatch_output_%j
#SBATCH --error=./sbatch_error_%j
source ~/.bash_profile
##########JOBSTART###########################

按照如上设定，编写文件运行job_somatic.pl这个文件。

第一次尝试：

 perl job_somatic.pl ./somatic_script/design.txt ./somatic_script/job_somatic.sh  32 hg19 /home/xxzhang/workplace/QBRC/geneome/hg19/hg19.fa /home/xxzhang/workplace/software/java/jdk1.7.0_80/bin/java 1 3 /home/xxzhang/workplace/QBRC/disambiguate_pipeline

Can’t exec “sbatch”: No such file or directory at job_somatic.pl line 54, line 3.

将第54行中的sbatch 修改为了bash之后，出现了新的报错：

batch accepts no parameters

尝试在example.sh中写一些内容。

perl job_somatic.pl example.txt job_somatic.sh 32 hg19 /home/xxzhang/workplace/QBRC/geneome/hg19/hg19.fa /home/xxzhang/workplace/software/java/jdk1.7.0_80/bin/java 3

傻傻的尝试，没有找对方法。后来经过多方的验证，我发现原来sbatch是slurm作业调度系统的作业提交的指令。而我们的服务器的作业调度系统是PBS，我应该根据我们的平台进行调整（所以跟着优秀的作者，还是能够学到很多有意思的东西的，学到新的东西真的很快乐）。

#!/bin/bash
# file: job_SNV.pbs
### set job name
#PBS -N example-job
### set output files
#PBS -o example.stdout
#PBS -e example.stderr
### set queue name
#PBS -q example-queue
### set number of nodes
#PBS -l nodes=2:ppn=4
# enter job's working directory
cd $PBS_O_WORKDIR
# get the number of processors
NP=`cat $PBS_NODEFILE | wc -l`
# run an example mpi4py job
mpirun -np $NP -machinefile $PBS_NODEFILE python example_mpi4py.py

保存文件名为job_SNV.pbs
而提交该文件的指令为，

qsub job_SNV.pbs

很粗糙的将job_somatic.pl中的sbatch修改为qsub之后，出错就改变了。说明我们的操作是有效的。

Use of uninitialized value $n in modulus (%) at job_somatic.pl line 30, line 1.
Illegal modulus zero at job_somatic.pl line 30, line 1.

剩下的就是，继续精细化的修改我们的job_somatic.sh指令。
另外一点，可能比较迷惑的是，为什么perl无法创建文件，somatic_calling_1.sh类如这样的文件。

下次写实验记录的时候，在开头的时候要清楚的总结下来，我们的这篇文章主要做了哪些方面的内容，这样的话，下次查找的时候比较容易查找。