实验记录 | 8/15-8/16 完成数据的合并和拆分

grep "AAACAGCCAATGAAGC" R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq | wc -l
67004
grep "AAACAGCCAATTGCGC" R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq | wc -l
48953

##################################################################################################

grep "GCCTCCATCCTCACAC" R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq | wc -l
597567
grep "TAAGTAGCACACCAAC" R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq | wc -l
343949

问题出现了！也许并不是我错了。
我的问题描述是这样的：
（1）在我的数据中数据量排名前几的数据，在师兄的list中找不到。
（2）在师兄的list中能够找到的，最大的，call出来的突变很少，且几乎都在chrM(按照学姐的说法，这个是核，mDNA的出现是一种干扰，没啥意义）

还是要再确认一下：

(base) [xxzhang@cu08 GEXst]$ grep "TCCGGAATCACACAGT" R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq |wc -l
152130
(base) [xxzhang@cu08 GEXst]$ grep "TGCTCCGTCTTAGCCC" R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq |wc -l
161069
(base) [xxzhang@cu08 GEXst]$ grep "CCGGTAGGTAAGCACC" R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq |wc -l
146639

另外拆分数据报错：

FileNotFoundError: [Errno 2] No such file or directory: ‘1_R2_chunk0261’
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File “scSplitter.py”, line 496, in
main()
File “scSplitter.py”, line 332, in main
pool.starmap(af.readBarcode, inputlist)
File “/home/xxzhang/miniconda3/lib/python3.6/multiprocessing/pool.py”, line 274, in starmap
return self._map_async(func, iterable, starmapstar, chunksize).get()
File “/home/xxzhang/miniconda3/lib/python3.6/multiprocessing/pool.py”, line 644, in get
raise self._value
FileNotFoundError: [Errno 2] No such file or directory: '1_R2_chunk0261’

这个错误，可咋解释呢？就是说找不到这个文件。但是，我实在不太明白chunk在这里是什么意思？为什么有的有chunk的R1、R2以及R5,有的，却没有呢？

#主要想验证一下，是不是像自己想的。因为整合之后会出现两个一模一样的序列，因此在后来拆分的时候报错。
(base) [xxzhang@cu08 combine_all]$ grep "@A00583:531:HW5FTDSXY:2:1101:1190:1016 1:N:0:GAGACGCA+ATGTTCAT" R21010305-3-wenku-S1-3-wenku-S1_all_R1.fastq
@A00583:531:HW5FTDSXY:2:1101:1190:1016 1:N:0:GAGACGCA+ATGTTCAT
#只有这一个#看来我的怀疑是错的！

#处理大文件一个感人的心得就是，要设计预实验，如果一开始就大数据的投资去做，一旦一步考虑错了，会有时间和金钱成本的巨大损耗。
#到现在都没有grep完成，就让他在这里跑着。
#其实现在遇到的主要的问题是：
#（1）即使是使用数据量最大的那几个细胞的fastq文件，得到的突变的数目依旧是非常的少。这就让我非常的迷茫。我猜想主要的原因是这样：
# * 我们使用的是rna的数据，rna在特定的细胞类群中表达。其本身的覆盖率是很大的问题。==>所以，我现在在考虑将第一次测的数据和第二次补测的数据整合，然后再拆成单个细胞，但是这块有400G,运行的特别慢。
# * 我们本身的10X，是一种非全长的测序，也就意味着范围又缩小了一半。
# * 我们使用的是核，而相较于核而言，线粒体的突变的累积会更多。
# * 我们使用的细胞是正常的神经组织的细胞，本身分裂的次数有限。因为突变主要是在分裂过程中累积下来的，故而连产生的这个环境都没有。
# * 我们使用的算法本身strelka，是用于寻找germline variant，而且是tumor-only的模式，是否需要更换，对于低vaf的variant的敏感性和特异性会更好？==>这一块需要文献调研？
# * 考虑到时间和效率的平衡，我们在做之前需要对我们要研究的数据进行初步的筛选？就是在实验设计这块下功夫？最终期待的结果是什么？

#总的来说，觉得自己对这个领域了解的还太少了。想多读些文献，看看别人都干了啥，进展到哪里？
#关于这个课题，我说不会放弃的。
#明天就跟老师谈论这些吧？看看怎么说？

---

今天，跟学姐吐槽之后，我发现了一件特别重要的事情。我觉得我自己对自己的要求要高一点。我明白了做事情应该要特别的严谨，不能够稀里糊涂的。做科研的基本素养，不能够弄虚作假，要大胆假设，小心求证。
我现在理一下，我现在还不够明确的地方：
（1）10X 数据的拆分问题，我现在还没有完全的搞定它。1000多行，看得我望而却步。心理上怕麻烦，但是实际上却是敷衍的态度，这并不好。
==>导致我现在对于这样的拆分的结果并不是特别的有自信，我甚至想要修改它的代码，但是我目前能力并不够。（我或许可以求助星宇帮忙）
（2）其次，我对于10X的数据的筛选的标准并不是特别的特别的自信。
==>因为10X的数据，由于技术层面的原因，而且我们又是正常的细胞，对于mutation calling的那个pipeline，究竟怎样的突变是germline究竟何者为somatic，不能够很好的区别？（我对于caller的基本原理望而生畏，是否适用于我们的数据？）
（3）对于SClineagar的参数的选择，目前只是使用他的默认参数在做，但是不知道需要修改些什么？
==>拒绝傻瓜式的被动的接受？

---

(base) [xxzhang@cu08 GEXst]$ head R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq
@A00583:531:HW5FTDSXY:2:1101:1190:1016 1:N:0:GAGACGCA+ATGTTCAT
GATTACTCATCACAGCTTGCACCAGGGTTTTTTTGTTTTTTTTTTTTTTTTTTTTTTTTTGTTGGGGTTGTTATGATTGACATCTGCAGCAGAGGGGTGACCATATCATACATATTTATATAAACTTTTTTTTTAAATATATATTTTTTT
+
F:FFF,F:,FFFFFF::FFFFFFF:,FF,,:,,F,FFFFFFFF:FF,F,,,F,,FF:F:F,FFFF,,,,FF,,,FF,,F,:,,,F:,,,:,,,,:F:F,F,,,,,,,F,F,F,,,F,,:,:FF,FF,,FFFF,,,F,,,:,,,FFFFFF,
@A00583:531:HW5FTDSXY:2:1101:1334:1016 1:N:0:GAGACGCA+ATGTTCAT
GCCTCCATCCTCACACACGGCTTTAGGCTTTTTTTTTTTTTTTTTTTTTTTTTGAGATACGGAGTCTGTCTCTGTTTCCCAGGCCGGAGTGCACTGGCCTGATCTCATCTCACTTCCACCTCTGCTTCCCTATTTCAAACAATTCTCCTT
+
FFF:F:FFF,FFFFFFFFFFFF,:FFFFFFFFFFFFFFFFFFFFFFFFFF,FFF:FF,FF,:FF,FFFFFFFFFFF,F,:FFF,F:FF:::F:FF,FF:F:FFFFF,,F:FF:F,,F,F:FF:,FF:FFF,,,,FF,:FFF,,:F:FFF,
@A00583:531:HW5FTDSXY:2:1101:1515:1016 1:N:0:GAGACGCA+ATGTTCAT
CGCATTTGTCCTAACTATGCTGTCAGCATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGTCCGCACACAACACCCAAACTCTCCCCAAAGATAAACCCCCCATTCAATAATACCTCTCTATTTTTCATTGTTTTATCAATATTTATTT
(base) [xxzhang@cu08 GEXst]$ cd ..
(base) [xxzhang@cu08 fastq]$ cd GEXnd/
(base) [xxzhang@cu08 GEXnd]$ head R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq
@A00265:638:H2HN7DSX2:4:1101:1063:1016 1:N:0:GAGACGCA+ATGTTCAT
ANGATGGAGCATGTCGTTCAATACCCATGGGGGGGGGGGGGGGCGGGTTTTTTTTTTCGCTCGGCCGCGGGCACCGGCTGCCCGCGCCCGTCCCGTCGGCCCCGCGGACGCCCCACTCTCCTGTCCGGCATCCCCCTGTTCAACTGCTGT
+
F#FFFFFFFFFFFFFFF:FFFFFFFFFFF::FF,FFFF,FFFF,FFF,,,F:F,,:FF,:,:,:,F,,F:F,,FFF::FF,:FFF,FF,,,,:FFF:FF,,F,F,,,,::,:,F,:,F,FF:,,F:F,F,F:F:FFFF,::::F,FF:::
@A00265:638:H2HN7DSX2:4:1101:1081:1016 1:N:0:GAGACGCA+ATGTTCAT
CNTGCGCGTCCCGAAGTCCGTGTACCCTCCAGCCCGGGCGGGGCGGGGTTTTTTTTCCCCGGACACGCGTCGGGCCCCGCGGGCCCCTCCGCGAAGAGCAGCGAGGGGGCGCCGCTAGGCACGTGGCGGGGCCGGTGGGTGGCTGGCCCC
+
F#FFFFFFFFFFFFFFFFFFFFFFFFFF,,,FF:F,FF,,FFFF,F:F:,F,,,,:FFF,,,,FF:,FF,,,F:,FFF,F,,::FFF,F:,,,,:F,F,,,,,,,FFF:::F,::F,FF:,F:,FF,,F:F:FF,,,:,,,:,,,:,:,,
@A00265:638:H2HN7DSX2:4:1101:1136:1016 1:N:0:GAGACGCA+ATGTTCAT
CNAGGGAGTCTAACCTTATCCCATTCTATGGTTTGTGGTTGGTGTTTTTTTTTTTACAAACGTTAATTTAAACACATATAATAAATTCCAACTAAACAAAAACCAGAAACAATAAATTTTTTTTCAAAAATTTTATTTAAAAATACAAAC

#查找在另一个数据中的title是否在同一个数据中出现
grep "@A00583:531:HW5FTDSXY:2:1101:1190:1016 1:N:0:GAGACGCA+ATGTTCAT" R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq
#没有结果

#查找同一个序列的sequence是否在另一个数据中出现 #我觉得既然是二次的重测，应该会找到一模一样的sequence的
grep 
"GATTACTCATCACAGCTTGCACCAGGGTTTTTTTGTTTTTTTTTTTTTTTTTTTTTTTTTGTTGGGGTTGTTATGATTGA               CATCTGCAGCAGAGGGGTGACCATATCATACATATTTATATAAACTTTTTTTTTAAATATATATTTTTTT"  
R21010305-3-wenku-S1-3-wenku-S1_combined_R1.fastq
#没有结果 #和我预想还真的不太一样

---

今天主要的任务，还是验证一下这个10X的数据是否能够很好的拆分？
刚刚拆了一下，觉得可以。最终还是投降了。现在选择宽松的拆第二次补测的数据。

---

我觉得会使用pbs之后，就觉得很方便了，特别的开心。算是这一周最开心的事情了。
但是刚刚qstat -n之后，发现别人好像比我跑的快。

mu01:
                                                                                  Req'd       Req'd       Elap
Job ID                  Username    Queue    Jobname          SessID  NDS   TSK   Memory      Time    S   Time
----------------------- ----------- -------- ---------------- ------ ----- ------ --------- --------- - ---------
10532.mu01              sytong      batch    jupyter          249270     1      4       --  7200:00:0 R 863:52:20
   cu08/25-28
10647.mu01              yqzhou      batch    rstudio          242092     1     32     200gb 7200:00:0 R 658:39:07
   cu03/0-31
10833.mu01              hywang      batch    rstudio          365009     1     16      50gb 7200:00:0 R 516:06:14
   cu06/0-15
11058.mu01              plli        batch    jupyter          342970     1     32       --  7200:00:0 R 295:01:31
   cu01/0-31
11145.mu01              wzhou       batch    STDIN            450522     1     32       --  7200:00:0 R 169:42:10
   cu07/0-31
11155.mu01              rqzhang     batch    jupyter.pbs      121712     1      8       --  7200:00:0 R 163:36:03
   cu08/3-10
11212.mu01              yqzhou      batch    STDIN             89722     1      1       --  7200:00:0 R  75:24:08
   cu08/0
11274.mu01              yuezhu      batch    chrom3diteration  45016     1      8       --  7200:00:0 R  56:07:40
   cu06/24-31
11281.mu01              whe         batch    jupyter           17244     1     25       --  240:00:00 R  49:45:14
   cu02/0-24
11284.mu01              yuezhu      batch    Jupyter          371613     1      8      50gb 7200:00:0 R  48:56:17
   cu08/12-19
11293.mu01              plli        fat      STDIN            151189     1     63       --  7200:00:0 R  30:48:56
   fat02/1-63
11306.mu01              plli        fat      STDIN            243373   --     --        --  7200:00:0 R  07:25:42
   fat02/0
11310.mu01              yuezhu      batch    chrom3diteration 260743     1      8       --  7200:00:0 R  04:43:27
   cu06/16-23
11312.mu01              yuezhu      batch    chrom3diteration 163311     1      8       --  7200:00:0 R  04:27:52
   cu08/1-2,11,20-24
11314.mu01              xxzhang     batch    split_10X_seq_2  181229     1      1       --  7200:00:0 R  00:19:49

我这里不太了解的有3个专业的术语。
（1）NDS：Resource_List.nodect = 1 #这个应该指的是使用的节点数
（2）TSK： Resource_List.nodes = 1:ppn=63 #这个应该指的是使用的多少核
（3）Queue：fat/batch的选择有何差别？

有的区别大了：

我下次可以学习一下。
我觉得这个实验室有很多我可以学习到的东西，尤其是向师兄师姐们学习。
我想我的进步是很快的，有一天我可以和他们一样的厉害。