对于所有数据类型，”文库拆分”涉及从一端或双端短序列中识别和移除细胞条形码序列 (cell barcode)和UMI。如果提前知道加入的细胞条形码，比如数据来自基于PCR板的方案，只需要找到条形码并与条形码库作比对，归类于与之最相似的那个就可以 (根据条形码的设计，一般允许最多1-2个错配)。这些数据通常在比对之前先做拆分，从而可以并行比对，提高效率。
我们可以利用很多公开的脚本，可以拆分任何plate-based的建库方案生成的数据，不管有没有UMI。
===================scRNASeqPipeline=========================
下载地址如下：https://github.com/tallulandrews/scRNASeqPipeline
perl 1_Flexible_UMI_Demultiplexing.pl read1.fq read2.fq b_structure index mismatch prefix
read1.fq : barcode/umi containing read
read2.fq : non-barcode containing read
b_structure : a single string of the format C##U# or U#C##
where C## is the cell-barcode and U# is the UMI.
e.g. C10U4 = a 10bp cell barcode followed by a 4bp UMI
index : file containg a single column of expected cell-barcodes.
if equal to \"UNKNOWN\" script will output read counts for each unique barcode.
mismatch : maximum number of permitted mismatches (recommend 2)
prefix : prefix for output fastq files.
perl 1_Flexible_UMI_Demultiplexing.pl read1.fq read2.fq C16U12 10cells_barcodes.txt 2 test
注：其实我在想这个index是指包含期望的cell barcode的序列文件，这个需要自己通过读原始数据去获得吗？还要用户自己提供？他们在读取输入的时候为啥不自己提取？
结果：对于index提供的没一个cell barcode生成一个文件，文件里面是包含这个cell barcode的测序reads。（为啥只有R2有结果，我R1里面明明除了16bp的cell barcode，12bp的UMI，还有150-16-12bp的reads，为什么不提取？我自己去运行cellranger的时候，发现它里面的bam文件也是如此，好像只比对了R2，没有处理R1）
=========================UMI tools====================================
下载地址：https://github.com/CGATOxford/UMI-tools
首先看下这个软件的功能：
UMI-tools功能模块
第一步：提取cell barcode白名单
whitelist 命令会从原始数据种提取去可能的cell barcode。通常情况下，10X的barcode长度为16nt，umi长度为10nt（但是我们的数据目前是16+12）；Drop-seq的barcode长度为12nt，umi长度为8nt。示例代码如下：
umi_tools whitelist -I pbmc_1k_v3_S1_L001_R1_001.fastq -p CCCCCCCCCCCCCCCCNNNNNNNNNNNN --plot-prefix=test --log2stderr > whitelist.txt
whitelist生成的结果文件包含四列：1、可接受的cell barcode；2、与可接受的barcode距离相差1的barcode；3、第一列barcode的umi数；4、第3列barcode的umi数
whitelist结果
同时也输出了真正的barcode的评估曲线，他们利用一个knee的概念，其实我也不懂这个。
barcode VS counts
然后利用找到的barcode进行提取和过滤相应的reads：
extract命令会从fastq文件中提取包含可接受barcode的reads，默认情况下extract命令会忽略umi的reads质量情况而不做处理。示例代码如下：
umi_tools extract --bc-pattern=CCCCCCCCCCCCCCCCNNNNNNNNNN --stdin hgmm_100_R1.fastq.gz --stdout hgmm_100_R1_extracted.fastq.gz --read2-in hgmm_100_R2.fastq.gz --read2-out=hgmm_100_R2_extracted.fastq.gz --whitelist=whitelist.txt
网友评论