实验记录 | 6/22 至 6/29

nohup perl somatic.pl NA NA ./example/example_dataset/sequencing/SRR7246238_1.fastq.gz ./example/example_dataset/sequencing/SRR7246238_2.fastq.gz 32 hg38 ./geneome/hg38/hg38.fasta /home/xxzhang/workplace/software/java/jdk1.7.0_80/bin/java ./output human 1 ./disambiguate_pipeline>DNA_pipeline_hg38.txt &

ERROR MESSAGE: Input files /home/xxzhang/workplace/QBRC/./geneome/hg38/hg38.fasta_resource/Mills_and_1000G_gold_standard.indels.hg38.vcf and reference have incompatible contigs.
ERROR /home/xxzhang/workplace/QBRC/./geneome/hg38/hg38.fasta_resource/Mills_and_1000G_gold_standard.indels.hg38.vcf contigs = [chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr20, chr21, chr22, chrX, chrY, chr1_KI270706v1_random, chr4_GL000008v2_random, chr14_GL000009v2_random, chr14_KI270726v1_random, chrUn_KI270742v1, chr1_KI270766v1_alt, chr7_KI270803v1_alt, chr8_KI270811v1_alt, chr8_KI270821v1_alt, chr14_KI270846v1_alt, chr15_KI270850v1_alt, chr17_KI270862v1_alt, chr17_KI270857v1_alt, chr22_KI270875v1_alt, chr22_KI270879v1_alt, chr2_KI270894v1_alt, chr17_KI270909v1_alt, chr19_KI270938v1_alt]
ERROR reference contigs = [chr10, chr10_GL383545v1_alt, chr10_GL383546v1_alt, chr10_KI270824v1_alt, chr10_KI270825v1_alt, chr11, chr11_GL383547v1_alt, chr11_JH159136v1_alt, chr11_JH159137v1_alt, chr11_KI270721v1_random, chr11_KI270826v1_alt, chr11_KI270827v1_alt, chr11_KI270829v1_alt, chr11_KI270830v1_alt, chr11_KI270831v1_alt, chr11_KI270832v1_alt, chr11_KI270902v1_alt, chr11_KI270903v1_alt, chr11_KI270927v1_alt, chr12, chr12_GL383549v1_alt, chr12_GL383550v2_alt, chr12_GL383551v1_alt, chr12_GL383552v1_alt, chr12_GL383553v2_alt, chr12_GL877875v1_alt, chr12_GL877876v1_alt, chr12_KI270833v1_alt, chr12_KI270834v1_alt, chr12_KI270835v1_alt, chr12_KI270836v1_alt, chr12_KI270837v1_alt, chr12_KI270904v1_alt, chr13, chr13_KI270838v1_alt, chr13_KI270839v1_alt, chr13_KI270840v1_alt, chr13_KI270841v1_alt, chr13_KI270842v1_alt, chr13_KI270843v1_alt, chr14, chr14_GL000009v2_random, chr14_GL000194v1_random, chr14_GL000225v1_random, chr14_KI270722v1_random, chr14_KI270723v1_random, chr14_KI270724v1_random, chr14_KI270725v1_random, chr14_KI270726v1_random, chr14_KI270844v1_alt, chr14_KI270845v1_alt, chr14_KI270846v1_alt, chr14_KI270847v1_alt, chr15, chr15_GL383554v1_alt, 等。

这个错误，我们之前也遇到过。现在我们对我们的vcf文档进行sortvcf的操作。

/home/xxzhang/workplace/software/java/jdk1.8.0_291/bin/java -jar /home/xxzhang/workplace/QBRC/somatic_script/picard.jar SortVcf --INPUT Mills_and_1000G_gold_standard.indels.hg38.vcf --OUTPUT Mills_and_1000G_gold_standard_sorted.indels.hg38.vcf -SEQUENCE_DICTIONARY /home/xxzhang/workplace/QBRC/geneome/hg38/hg38.dict

报出错误：

SAM dictionaries are not the same: SAMSequenceRecord(name=chr1,length=248956422,dict_index=0,assembly=20,alternate_names=[]) was found when SAMSequenceRecord(name=chr10,length=133797422,dict_index=0,assembly=null,alternate_names=[]) was expected.

最后找出罪魁祸首。我们的参考数据集本身就有问题。

chr10
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

为什么chr10是最开始的呢？这明显不对。
于是我们换用从服务器上下载到的hg38的序列看一下。
参考基因组的序列：

chr1 248956422 112 100 101
chr2 242193529 251446211 100 101
chr3 198295559 496061788 100 101
chr4 190214555 696340415 100 101
chr5 181538259 888457250 100 101
chr6 170805979 1071811004 100 101
chr7 159345973 1244325155 100 101
chr8 145138636 1405264700 100 101
chr9 138394717 1551854835 100 101
chr10 133797422 1691633613 100 101
chr11 135086622 1826769123 100 101
chr12 133275309 1963206725 100 101
chr13 114364328 2097814901 100 101
chr14 107043718 2213322999 100 101
chr15 101991189 2321437268 100 101
chr16 90338345 2424448481 100 101
chr17 83257441 2515690322 100 101
chr18 80373285 2599780450 100 101
chr19 58617616 2680957593 100 101
chr20 64444167 2740161498 100 101
chr21 46709983 2805250232 100 101
chr22 50818468 2852427440 100 101
chrX 156040895 2903754205 100 101
chrY 57227415 3061355656 100 101
chrM 16569 3119155468 100 101
chr1_KI270706v1_random 175055 3119172340 100 101
chr2_KI270715v1_random 161471 3119923517 100 101
chr2_KI270716v1_random 153799 3120086740 100 101
chr3_GL000221v1_random 155397 3120242214 100 101
chr4_GL000008v2_random 209709 3120399302 100 101
chr5_GL000208v1_random 92689 3120611245 100 101
chr9_KI270717v1_random 40062 3120704997 100 101
chr9_KI270718v1_random 38054 3120745596 100 101
chr9_KI270719v1_random 176845 3120784168 100 101
chr9_KI270720v1_random 39050 3120962918 100 101
chr11_KI270721v1_random 100316 3121002498 100 101
chr14_GL000009v2_random 201709 3121103957 100 101
chr15_KI270727v1_random 448248 3122208617 100 101
chr16_KI270728v1_random 1872759 3122661488 100 101
chr17_GL000205v2_random 185591 3124553114 100 101
chrY_KI270740v1_random 37240 3126188990 100 101
chrUn_KI270302v1 2274 3126226720 100 101
chrUn_KI270304v1 2165 3126229134 100 101
chrUn_KI270303v1 1942 3126231438 100 101

我们在注释文件vcf中的序列顺序：

##contig=<ID=chr1,length=248956422,assembly=20>
##contig=<ID=chr2,length=242193529,assembly=20>
##contig=<ID=chr3,length=198295559,assembly=20>
##contig=<ID=chr4,length=190214555,assembly=20>
##contig=<ID=chr5,length=181538259,assembly=20>
##contig=<ID=chr6,length=170805979,assembly=20>
##contig=<ID=chr7,length=159345973,assembly=20>
##contig=<ID=chr8,length=145138636,assembly=20>
##contig=<ID=chr9,length=138394717,assembly=20>
##contig=<ID=chr10,length=133797422,assembly=20>
##contig=<ID=chr11,length=135086622,assembly=20>
##contig=<ID=chr12,length=133275309,assembly=20>
##contig=<ID=chr13,length=114364328,assembly=20>
##contig=<ID=chr14,length=107043718,assembly=20>
##contig=<ID=chr15,length=101991189,assembly=20>
##contig=<ID=chr16,length=90338345,assembly=20>
##contig=<ID=chr17,length=83257441,assembly=20>
##contig=<ID=chr18,length=80373285,assembly=20>
##contig=<ID=chr19,length=58617616,assembly=20>
##contig=<ID=chr20,length=64444167,assembly=20>
##contig=<ID=chr21,length=46709983,assembly=20>
##contig=<ID=chr22,length=50818468,assembly=20>
##contig=<ID=chrX,length=156040895,assembly=20>
##contig=<ID=chrY,length=57227415,assembly=20>
##contig=<ID=chrM,length=16569,assembly=20>
##contig=<ID=chr1_KI270706v1_random,length=175055,assembly=20>
##contig=<ID=chr1_KI270707v1_random,length=32032,assembly=20>
##contig=<ID=chr1_KI270708v1_random,length=127682,assembly=20>
##contig=<ID=chr1_KI270709v1_random,length=66860,assembly=20>
##contig=<ID=chr1_KI270710v1_random,length=40176,assembly=20>
##contig=<ID=chr1_KI270711v1_random,length=42210,assembly=20>
##contig=<ID=chr1_KI270712v1_random,length=176043,assembly=20>
##contig=<ID=chr1_KI270713v1_random,length=40745,assembly=20>
##contig=<ID=chr1_KI270714v1_random,length=41717,assembly=20>
##contig=<ID=chr2_KI270715v1_random,length=161471,assembly=20>
##contig=<ID=chr2_KI270716v1_random,length=153799,assembly=20>
##contig=<ID=chr3_GL000221v1_random,length=155397,assembly=20>
##contig=<ID=chr4_GL000008v2_random,length=209709,assembly=20>
##contig=<ID=chr5_GL000208v1_random,length=92689,assembly=20>
##contig=<ID=chr9_KI270717v1_random,length=40062,assembly=20>
##contig=<ID=chr9_KI270718v1_random,length=38054,assembly=20>
##contig=<ID=chr9_KI270719v1_random,length=176845,assembly=20>
##contig=<ID=chr9_KI270720v1_random,length=39050,assembly=20>
##contig=<ID=chr11_KI270721v1_random,length=100316,assembly=20>
##contig=<ID=chr14_GL000009v2_random,length=201709,assembly=20>
##contig=<ID=chr14_GL000225v1_random,length=211173,assembly=20>
##contig=<ID=chr14_KI270722v1_random,length=194050,assembly=20>

总体而言，这两者的顺序是一致的。我们可以试一试！
后来，成功。
另外，除了之前，我们建立的索引文件以及hg38_resource的文件之外，我们还需要对annovear中的注释文件进行更新。
更新完成之后，我们再进行分析。

---

在这个间隙，我们调查一下，我们复现过程中需要用到的数据集都有哪些？

（1）The SMART-Seq datasets of 969 CD45+ cells

PMID: 29884728 (Wang et al., 2018a)
EGA: EGAS00001000926

数据检索链接：https://ega-archive.org/search-results.php?query=EGAS00001000926

在引用的原始的文献中，对于数据的来源是这样描述的：

Data are provided for UTSW patients specifically consenting to placement of their raw genomic data in a protected publicly accessible database. Sequencing data from consented patients are deposited in EGA, with accession numbers EGAS00001002786 and EGAS00001000926.

这个文献的作者还是那个王涛。

tao.wang@utsouthwestern.edu

已发送邮件询问，静待结果（21：01）。得到回答，这个回答并不能提供有效的信息（21：10）。

这个数据的作者的联系方式是：

gdac-rcc-d@gene.com

需要弄清楚作者使用的数据是哪一张？
数据连接：https://ega-archive.org/studies/EGAS00001000926

（2）The SMART-Seq datasets of CML cells

PMID: 28504724 (Giustacchini et al., 2017)
GEO： GSE76312
SRA： SRP067759

数据检索链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76312
原始数据链接：https://www.ncbi.nlm.nih.gov/sra?term=SRP067759

这个数据的分类比较复杂。主要使用的是CML的一个样本绘制的那张热图。===>近期目标。
不过更加让人无法抉择的是，参考比对所使用的normal样本，作者是如何选择的？

我们就以CML15的诊断的patient的样本为例，来绘制热图。

===>经过与作者的讨论，我们得到了作者使用的样本的具体信息。
首先，作者使用K562样本绘制VAF的热图，然后从中抽取出突变ChrX 12975141 T->A、Chr6 132814747 A–>G、Chr6 35469489 G–>GA,绘制barplot的图。

其次，作者使用K562、CML1266、CML656、OX1931这四个样本（患者）的数据，绘制下面的boxplot的图。

===>明确了这些信息之后，我更加清晰的知道自己该做什么了。

我先看一下这四个patient是什么样的数据特征：

Authenticated K562 and mycoplasma-negative (chronic myeloid leukemia, human cell line) cells were obtained from the American Type Culture Collection (ATCC) and grown in Iscove’s modified Dulbecco’s media (IMDM), 10% fetal bovine serum (FBS).
K562细胞是第一个人类髓性白血病人工培养的细胞。红白血病K562细胞的类型,是源自一个53岁的女性慢性髓性白血病爆发期病人的淋巴母细胞。用于研究肿瘤和白血病治疗、药物靶标等领域。（所以，我才知道了这是怎么回事。这个细胞系是一个医学上常用的白血病细胞系。）

CML1266
bcr-abl status: negative
tumor stage: pre_blast_crisis
phenotype: Lin-CD34+CD38-

CML656
bcr-abl status: negative
tumor stage: blast_crisis
phenotype: Lin-CD34+CD38-

OX1931
bcr-abl status: positive
tumor stage: pre_blast_crisis
phenotype: Lin-CD34+CD38-

我现在已经使用同样的命令，将数据下载完成。下载完成之后，我们需要将数据从sra==>fastq.gz。
根据网站https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=toolkit_doc&f=fastq-dump
我们可能需要的参数有：

–split-files (我们的数据都是单端的，所以只是输出一条数据)
–outdir（输出文件夹）
–gzip （输出结果压缩）

所以，我们现在需要做的是编写批处理的指令。
先解决单个文件的问题。

fastq-dump --split-files --gzip ./SRR3052037/SRR3052037.sra --outdir ./SRR3052037/

现在的任务，就是想把其变为批处理的指令。
我去学习一下bash文件如何书写。

nohup fastq-dump --split-files --gzip /home/xxzhang/workplace/QBRC/data/CML/K562/SRR*/SRR*.sra &

按照这种方式，我们书写批处理的代码为：

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

程序正常的运行中，我还真的挺喜欢通过代码来处理问题的。安身立命，安身立命，我还挺喜欢通过代码来批量的，高效的，快速的处理问题的。既然有这么一个得力的帮手，我们为什么不好好的运用呢？这更加坚定了我要坚持学代码的决心。要把这个工具很棒的用起来。

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

（3）The 10X Genomics single cell sequencing

数据链接：https://www.10xgenomics.com/resources/datasets/10X
Genomics: Vdj_v1_hs_pbmc2_t; pbmc4k

我似乎找到了一直想要找到的，pumc4k的数据。现在正在下载中，
pbmc4k： https://cf.10xgenomics.com/samples/cell-exp/1.3.0/pbmc4k/pbmc4k_web_summary.html
https://www.10xgenomics.com/resources/datasets/4-k-pbm-cs-from-a-healthy-donor-2-standard

Vdj_v1_hs_pbmc2_t： https://cf.10xgenomics.com/samples/cell-vdj/3.0.0/vdj_v1_hs_pbmc2_t/vdj_v1_hs_pbmc2_t_web_summary.html
https://www.10xgenomics.com/resources/datasets/pbm-cs-of-a-healthy-donor-tcr-enrichment-from-amplified-c-dna-1-standard
找到数据之后，下载。

wget https://cf.10xgenomics.com/samples/cell-vdj/3.0.0/vdj_v1_hs_pbmc2_t/vdj_v1_hs_pbmc2_t_fastqs.tar
wget https://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/1.3.0/pbmc4k/pbmc4k_fastqs.tar

我们来看一下，这两套数据的详细信息。

下载完成之后，需要使用作者再文献中推荐的方法，将其拆分。
我们来看一下原文。

The 10x Genomics single-cell sequencing platform outputs the raw sequencing reads of all cells in one Fastq file. We developed ScSplitter to split the reads in original Fastq file by their cell of origin using their cell barcode sequences.

参考链接：https://github.com/zzhu33/scSplitter

wget https://github.com/zzhu33/scSplitter/raw/master/scSplitter_v0.9.5.zip
unzip  scSplitter_v0.9.5.zip 
cd scSplitter/

al_functions_m95.py exampleInputNames.txt scSplitter.py

tar -xvf pbmc4k_fastqs.tar 
tar -xvf vdj_v1_hs_pbmc2_t_fastqs.tar 

解压之后，文件夹中得到子文件。
程序如果要继续运行，需要STAR建立索引。不过目前我好想还没有hg38.gtf文件，因此，我需要先去服务器或者其他地方看一下。
没有。那么gtf文件，我们明天下载并处理完成。

/home/xxzhang/workplace/software/STAR/bin/Linux_x86_64/STAR --runThreadN 32 --runMode genomeGenerate --genomeDir /home/xxzhang/workplace/QBRC/geneome/hg38/STAR --genomeFastaFiles hg38.fasta --sjdbGTFfile hg38.gtf --sjdbOverhang 149

我差点忘记我今天的主要任务是什么了！（早上开组会开了超级久，一上午就过去了）
现在是下午（也快过完了），我要去找hg38的gtf文档。我查找一下我的实验记录，hg19的实验文档，我是如何寻找的？

 wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/genes/hg38.ensGene.gtf.gz
 gunzip hg38.ensGene.gtf.gz
 sort -k1,1V -k3,3n -k4,4n hg38.gtf >hg38_sorted.gtf
 grep -o -E "^\w+([-+.]\w+)*" hg38_sorted.gtf|uniq 
 mv hg38_sorted.gtf hg38.gtf

这样操作之后，就是对其进行构建索引。在索引的构建中。

索引构建完成之后，我们将其进行比对。

python3 scSplitter.py --ig True --f pbmc4kInputNames.txt  --i   /home/xxzhang/workplace/QBRC/data/10X/pbmc4k   --r  /home/xxzhang/workplace/QBRC/data/10X/pbmc4k_cells   --ind    /home/xxzhang/workplace/QBRC/geneome/hg38/STAR

报出错误：

ModuleNotFoundError: No module named ‘pandas’

pip3 install psutil
pip3 install pandas

出现错误。

ERROR -2:invalid number of input files: 1

那么，面对这种情况，该怎么办呢？
查看了文件的源代码，发现报错的环节是这样提示的：

if (n0 > 3) or (n0 < 2):
    print(f'\n\nERROR -2:invalid number of input files: {n0}\n\n')
    return -2

分析原因是，.txt文件必须是以tab键隔开的。于是，如上修正，得以顺利的运行。

出现了错误2：

FileNotFoundError: [Errno 2] No such file or directory: ‘lane_1Aligned.out.sam’

终于找到了出错的地方：

/bin/sh: STAR: command not found

虽然到了这一步，但是我不知道如何修改它的原始代码。问题的症结就在于，无法告知程序其STAR所在的位置。
星宇说，是不是需要提前配置环境变量！
我将STAR的路径配置进去，看看最后的结果，报不报错。

另外我想确认的是，这个数据的Chemistry的version是不是V2？
pbmc4k的Chemistry的version是：

经过图片显示是V3 Chemistry。
pbmc2的Chemistry的version是：

经过图片显示是V1。

那么这个时候就出现了另一个问题，Chemistry的版本对数据的处理结果影响大吗？
而为什么提出这个问题，是因为我们的代码中，有对于这个version的判断。

parser.add_argument('--chver', type=int, default=2,
                    help='10x chemistry version; options: 1 (v1 chemistry paired), 2 (v2 chemistry); default: 2')

parser.add_argument('--tsol', type=int, default=13,
                    help='length of TSO, for older 10x Genomics v1 chemistry paired end sequencing; default = 13 bp')
if chemVer == 1:   
    versionStr = '10x v1 chemistry paired'
if chemVer == 1:
    isV1chem = True
else:
    isV1chem = False
if chemVer == 1:
    totalReads *= 2
if isV1chem:
    af.alignFastq(read1str, f'{pathS}/lane_{i + 1}', refind, nprocess, read2str)
else:
    af.alignFastq(read2str, f'{pathS}/lane_{i + 1}', refind, nprocess)

该代码，对于此处有相当大的依赖性。

另外，又出错了。

IndexError: list index out of range
selecting reads based on alignment…
multiprocessing.pool.RemoteTraceback:
“”"
Traceback (most recent call last):
File “/home/xxzhang/miniconda3/lib/python3.6/multiprocessing/pool.py”, line 119, in worker
result = (True, func(*args, **kwds))
File “/home/xxzhang/miniconda3/lib/python3.6/multiprocessing/pool.py”, line 47, in starmapstar
return list(itertools.starmap(args[0], args[1]))
File “/home/xxzhang/workplace/software/scSplitter/al_functions_m95.py”, line 407, in readSelection
rname = list(zip(*saml))[2]
IndexError: list index out of range
“”"

我觉得出错的主要原因是，默认值的问题，这个值的选择应该基于我们的数据的基本情况。

parser.add_argument(’–cb’, type=int, default=16,
help=‘length of cell barcode; default = 16 bp’)
parser.add_argument(’–sl’, type=int, default=8,
help=‘length of sample index; default = 8 bp’)
parser.add_argument(’–ul’, type=int, default=10,
help=‘length of UMI; default = 10 bp’)

这里说的，list index out of range,指的就是，index的默认值超出了边界。
我觉得之所以会出错，可能是chemistry的版本不同，而此处我们使用的是默认的index的8bp,而可能对于我们自己的数据而言，v3的情况会有所不同。
所以，对于这里的length of cell barcode/length of sample index/length of UMI这些概念，我们都需要认真的明确一下。

这张图显示的是，10X的chemistry v1的数据的主要的建库的信息。
而chemistry v2和v3的建库信息尚未在官网中查找到。
我现在想用scSplitter来对pbmc2的数据（v1 Chemistry）进行拆分。

python3 scSplitter.py --ig True --f  pbmc2kInputNames.txt  --i   /home/xxzhang/workplace/QBRC/data/10X/vdj_v1_hs_pbmc2_t   --r  /home/xxzhang/workplace/QBRC/data/10X/pbmc2_v1_cells   --ind    /home/xxzhang/workplace/QBRC/geneome/hg38/STAR

pbmc2kInputNames.txt文件的主要内容：

vdj_v1_hs_pbmc2_t_S1_L006_I1_001.fastq.gz vdj_v1_hs_pbmc2_t_S1_L006_R1_001.fastq.gz vdj_v1_hs_pbmc2_t_S1_L006_R2_001.fastq.gz
vdj_v1_hs_pbmc2_t_S1_L007_I1_001.fastq.gz vdj_v1_hs_pbmc2_t_S1_L007_R1_001.fastq.gz vdj_v1_hs_pbmc2_t_S1_L007_R2_001.fastq.gz

将指令切换到/home/xxzhang/workplace/software/scSplitter的路径下。

python3 scSplitter.py --ig  True --f  /home/xxzhang/workplace/QBRC/data/10X/pbmc2kInputNames.txt  --i  /home/xxzhang/workplace/QBRC/data/10X/vdj_v1_hs_pbmc2_t/  --r /home/xxzhang/workplace/QBRC/data/10X/pbmc2_v1_cells  --ind  /home/xxzhang/workplace/QBRC/geneome/hg38/STAR --chver 1 --gz TRUE

遇到错误：

FATAL ERROR in reads input: quality string length is not equal to sequence length.

根据网上的介绍（参考链接：https://groups.google.com/g/rna-star/c/kUI6KgkTJPg），大概率是fastq文件的问题。

each read is represented with 4 lines in the FASTQ files@readID
sequence
quality scorese.g.:@D00102:CBR8BANXX171122:CBR8BANXX:1:1101:10003:18149 1:N:0:TAAGGCGACTCTCTAT
GCTCGTAGGAGCGTATCATCAGCTGGGTGCCGTTCTTCCTGTCTCTTATAC
/B<</BFFFFFFFFBBFFFFFFFFFFFFFFFFBFFFFFFFFFBBFFBFFFF

The length of the quality string should be equal to the sequence length.就是说，在fastq文件中，有序列中每一个碱基的序列信息，也有序列中每一个碱基的测序质量。因此，两者的长度必须一致，否则如何一一对应呢？因此，在这个层面上，我要去看一下我的错误信息中，使用的fastq序列是什么？

–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

我首先怀疑是不是我的输入文件的问题，文件之间存在对应关系。
将数据调整了一下，R1，R2，I1的顺序。
重新参照之前的命令，书写代码：

python3 scSplitter.py --ig True --f  pbmc2kInputNames.txt  --i   /home/xxzhang/workplace/QBRC/data/10X/vdj_v1_hs_pbmc2_t   --r  /home/xxzhang/workplace/QBRC/data/10X/pbmc2_v1_cells   --ind    /home/xxzhang/workplace/QBRC/geneome/hg38/STAR  --chver 1  --gz TRUE

出现错误：

ERROR -2:invalid number of input files: 4

后来经过错误的排查，我们发现主要的原因在于，我的pbmc2kInputNames.txt 的文件中，有“看不见”的空格，因此，计算机比我们敏锐的多，报出错误。我排查的主要方式是修改 scSplitter.py文件，让其输出row值，通过输出的row值，我们发现了多余的空格。

the rows line is : [‘vdj_v1_hs_pbmc2_t_S1_L006_R1_001.fastq.gz’, ‘vdj_v1_hs_pbmc2_t_S1_L006_R2_001.fastq.gz’, ‘vdj_v1_hs_pbmc2_t_S1_L006_I1_001.fastq.gz’, ‘’]

因此将这个空格删去即可。
另外，程序接着运行的过程中，还报了其他的错误。我们在另一篇博客中继续去解决这个问题。

（4）The CTCL ECCITE-seq dataset

PMID: 31011186(Mimitou et al., 2019)
GEO: GSE126310

参考链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126310

这个数据是存放在dbGaP数据库中的。关于这块数据的使用，我之前有问过这个作者。这套数据集中，关于CTCL的样本有五个。

GSM3596100 CTCL-cDNA （需要申请）
GSM3596101 CTCL-ADT
GSM3596102 CTCL-HTO
GSM3596103 CTCL-TCRab（需要申请）
GSM3596104 CTCL-TCRgd（需要申请）

这五个样本中，

CTCL-cDNA is used for calling variants and inferring missing data points.
CTCL-TCRab is used for lineages validation
being submitted to dbGaP due to patient privacy concerns

需要去dbGaP数据库上申请。

（5）The Liver mCelSeq2 dataset

PMID: 31292543(Aizarani et al., 2019)
GEO：GSE84498
SRP：SRP078795

数据链接：https://www.ncbi.nlm.nih.gov/sra?term=SRP078795
这个数据好像就是小鼠的正常的样本，非tumor。这种情况下，应该如何选择文件呢？
我们先尝试批量的下载sra数据文件，然后再将其转化为fastq.gz的文件。

wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-centos_linux64.tar.gz
tar xzvf sratoolkit.current-centos_linux64.tar.gz
rm -f sratoolkit.current-centos_linux64.tar.gz
mv sratoolkit.2.11.0-centos_linux64/ sratoolkit
cd sratoolkit/
cd bin/
vi ~/.bashrc
source ~/.bashrc
vdb-config -i
prefetch --option-file /home/xxzhang/workplace/QBRC/data/mCel/SRR_Acc_List.txt --output-directory /home/xxzhang/workplace/QBRC/data/m                                                                                                             Cel/

其中SRR_Acc_List.txt的内容如下：

SRR3928573
SRR3928574
SRR3928575
SRR3928576
SRR3928577
SRR3928578
SRR3928579
SRR3928580
SRR3928581
SRR3928582
SRR3928583
SRR3928584
SRR3928585
SRR3928586
SRR3928587
SRR3928588
SRR3928589
SRR3928590
SRR3928591
SRR3928592
SRR3928593
SRR3928594
SRR3928595
SRR3928596
SRR3928597
SRR3928598

至于这个.txt文件的内容，我们是在SRA数据库中，通过右上角Send to==>File==>Format==>Runinfo==>Create File的步骤得到。即可，得到具体的.txt文件。

https://www.spatialresearch.org/resources-published-datasets/doi-10-1126science-aaf2403

数据也正在下载中。

（8）The scATAC-Seq datasets

PMID: 30827679 (Ludwig et al., 2019)
GEO: GSE115218

原来我以为没有写SRP的ID，后来才发现这个文件的原始数据还是可以获取的。虽然有2000+个样本，但是单个其实并不大（~10Mb）。
这个数据集中的类别还是挺多的，有bulk的数据，也有sc的数据，也有mtDNA的数据。而我们需要的就是，scATAC-seq的数据，所以需要筛查一下。

首先从SRP界面下载注释信息。
参考链接：https://www.ncbi.nlm.nih.gov/sra?linkname=bioproject_sra_all&from_uid=474183
其次到GEO accession的界面下载注释信息。
参考链接：https://www.ncbi.nlm.nih.gov/geo/browse/?view=samples&series=115218
把两个list整理完成了。接下来就是要处理怎样批量的运行产生结果了。

nohup ./prefetch --option-file /home/xxzhang/workplace/QBRC/data/scATAC/scATAC_seq_TF1.txt --output-directory /home/xxzhang/workplace/QBRC/data/scATAC/ &
nohup ./prefetch --option-file /home/xxzhang/workplace/QBRC/data/scATAC/scATAC_seq_CC100.txt --output-directory /home/xxzhang/workplace/QBRC/data/scATAC/ &