UMI-tools | Cellranger手动分析流程 | 超短reads（primer、barcode、UMI、index等）比对方法 | gRNA | 单细胞 | Perturb-seq
2023年03月31日
cellranger pipeline基本是对的，但想要更原始的数据还是得自己手动count。
这两批perturb-seq的问题在于，Plasmid的设计，无法区分原始的和编辑后的Plasmid，长度一样，于是做cell sorting的时候一大堆dummy cells就被筛选出来了，真正有gRNA的就很少。
以下是手动counting的代码，建议用linux命令，python很慢很慢很慢。
提取fastq里的序列
zcat R2_extracted.fastq.gz | sed -n '1~4p' | cut -f2 -d'_' > CB zcat 1_pBA904_Std_Perturb_seq_327809w_CH9.fastq.gz | sed -n '2~4p' > raw_gRNA
grep可以提取line
grep --line-buffered AAAAAGGGGAGAGAGG whitelist.txt grep --line-number AAAAAGGGGAGAGAGG whitelist.txt
批量运行
cat /home/zz950/tmpData/new_seq/crispr_reference_feature_pBA904.csv | while IFS="," read id name read pattern sequence type; do if [ $sequence = sequence ]; then continue fi echo $sequence grep --line-number $sequence raw_gRNA | cut -f1 -d':' > $id.grep done
grep --ignore-case --line-number ttggaccaggatgggcaccacccgtttaagagc raw_gRNA | cut -f1 -d':' > empty.grep
非常慢的python代码也贴出来吧：
import pysam
import pandas as pd
df = pd.DataFrame()
filename = "R2_extracted.fastq.gz"
fh = pysam.FastxFile(filename)
# with pysam.FastxFile(filename) as fh:
i = 0
for entry in fh:
i = i + 1
# print(entry.name)
# print(entry.sequence)
# print(entry.comment)
# print(entry.quality)
CB = entry.name.split("_")[1]
UMI = entry.name.split("_")[2]
gRNA = entry.sequence
# df = df.append({'CB': CB, 'UMI': UMI, 'gRNA_raw': gRNA}, ignore_index=True)
tmp_df = pd.DataFrame([[CB, UMI, gRNA]], columns=['CB', 'UMI', 'gRNA_raw'])
df = pd.concat([df, tmp_df])
if i % 1000 == 0:
print(i) # 66124793
# break
fh.close()
df.to_csv("CB_gRNA_raw.csv")
2023年03月28日
新一批perturb-seq出来了，cellranger的分析结果很离谱，暂时不知道是为什么，cell barcode和gRNA都可以通过grep找到，最终却什么都没有。
zcat CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R1_001.fastq.gz | grep AAAGAACAGGTCATAA
于是就对CRISPR lib做自己的处理，UMI-tools提取barcode，然后自己做gRNA count矩阵。
最好对R1做一个trim：
seqtk trimfq -b 14 -e 108 raw/CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R1_001.fastq.gz > CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R1_001.fastq seqtk trimfq -b 24 -e 90 raw/CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R2_001.fastq.gz > CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R2_001.fastq
# Step 1: Identify correct cell barcodes
umi_tools whitelist --stdin CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R1_001.fastq.gz \
--bc-pattern=CCCCCCCCCCCCCCCCNNNNNNNNNN \
# --set-cell-number=100 \
--log2stderr > whitelist.txt
# Step 2: Extract barcdoes and UMIs and add to read names
umi_tools extract --bc-pattern=CCCCCCCCCCCCCCCCNNNNNNNNNN \
--stdin CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R1_001.fastq.gz \
--stdout R1_extracted.fastq.gz \
--read2-in CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R2_001.fastq.gz \
--read2-out=R2_extracted.fastq.gz \
--whitelist=whitelist.txt
参考：
2022年08月24日
上一次提需求是在2017年，在华大内部参与单细胞流程搭建，那时候什么都需要自己做，比如提取barcode和index，现在cellranger流程里都自带了。
现在则是在单细胞Perturb-seq分析时遇到了问题，提供的gRNA比对不佳，只有10%不到的reads有gRNA，这离理想的60%以上差太远了，必须得找原因了。
这里我想将我们的gRNA比对到fastq上，看到底是什么原因，我猜测是gRNA左右可能被截断了，所以你拿全长去比对肯定有问题，而且cellranger只允许一个mismatch，而且无参数可调。
cellranger其实是一个非常封闭的软件，非常呆板，碰到问题时很多个性化的功能都没有，好处是极大的加速了分析的步伐，缺点就是分析出来的结果质量不佳，会导致市场上的劣质结果大量堆积。
perturb-seq是一个非常新的领域，有些坑你不得不填，比如这个最核心的gRNA的比对问题，显然cellranger的提取是有问题的！
但不要随便深入一个算法开发，坑容易开，但很难填，一点要学会利用可靠的现成工具，通过整合它们来解决你遇到的问题。
UMI-tools - https://umi-tools.readthedocs.io/en/latest/index.html
Alevin - https://salmon.readthedocs.io/en/latest/alevin.html
终于找到使用方法了，用UMI-tools把R1的fastq的cell barcode和UMI提取到read name上，然后比对找出R2最佳匹配的gRNA，就这么简单。
-
https://umi-tools.readthedocs.io/en/latest/regex.html
- https://umi-tools.readthedocs.io/en/latest/QUICK_START.html
- https://divingintogeneticsandgenomics.rbind.io/post/understand-10x-scrnaseq-and-scatac-fastqs/
- Analysing 10X Single Cell RNASeq Data
10x的技术问题：
umi_tools extract [OPTIONS] -p PATTERN [-I IN_FASTQ[.gz]] [-S OUT_FASTQ[.gz]] --read2-in=IN2_FASTQ[.gz] --read2-out=OUT2_FASTQ[.gz] umi_tools extract --whitelist=/home/lizhixin/softwares/cellranger-7.0.1/lib/python/cellranger/barcodes/3M-february-2018.txt --stdin=HT29_P1_CRISPR_CKDL220019395-1A_H7MJYDSX5_S1_L003_R1_001.fastq.gz --bc-pattern=CCCCCCCCCCCCCCCCNNNNNNNNNN --log=processed.log --stdout processed.R1.fastq.gz --read2-in=HT29_P1_CRISPR_CKDL220019395-1A_H7MJYDSX5_S1_L003_R2_001.fastq.gz --read2-out=processed.R2.fastq.gz
bowtie-build -f feature_ref.fasta ref/gRNA
bowtie [options]* <ebwt> {-1 <m1> -2 <m2> | --12 <r> | <s>} [<hit>]
bowtie ./ref/gRNA -q <( gunzip < /home/lizhixin/project/Data_center/analysis/Perturb_seq_Sethi/rawData/fastqs/HT29_P1_CRISPR/processed.R2.fastq.gz ) --seedmms 2 --seedlen 18 --trim5 25 --trim3 70 > gRNA.sam
--tryhard
bowtie似乎没有好结果，我决定用R来做。
顺便画一下gRNA的motif，是否有很强的变异。
zcat processed.R2.fastq.gz | grep -F -B 1 -A 2 'ACCGACTTCCTCCGC|CAACACTGTCAATGT|ACCTGCTACGACAGT|CCGGGCAACGTTACG|AAGAGTCCATCAAAC|TCAATGCCATCCCGC|ACCAGGACAGGAACA|GGCTGCTGAACTTCC|ACAGAGGGCAAAAAA|ACTTCTGAGGCTCCA|GTTATGCAAGGTCCC|TGTGCTGTGATTGGCG|TCAGCGTCCGACCCAG|GATGCGCTCACCTCCC|TCTTGGCCATGGTGCC|GCGCCTGCGCACTGAG|GAGACGGAGGAAGGTC|GACTCCGGCCGGCTCT|CCGCGTGCGTCCCGGGTTAC|CCATCTGTCCTAATTCGGAT' > tmp.fastq conda install -c tsnyder tre zcat HT29_P1_CRISPR_CKDL220019395-1A_H7MJYDSX5_S1_L003_R2_001.fastq.gz | agrep -1 -n 'TCACCGACTTCCTCCGCCG' | wc -l
用R搞了一下午，发现比对并没有什么很好的结果，mismatch和indel的概率并不高，可能就是data本身有问题。
序列比对这个东西太细节太技术了，需要做大量测试，自己玩玩还行，真要可靠性那还是用现成的比对软件吧，否则自己的结果自己都不信。
这个才是正确的grep代码
zcat processed.R2.fastq.gz | grep -B 1 -A 2 'ACCGACTTCCTCCGC\|CAACACTGTCAATGT\|ACCTGCTACGACAGT\|CCGGGCAACGTTACG\|AAGAGTCCATCAAAC\|TCAATGCCATCCCGC\|ACCAGGACAGGAACA\|GGCTGCTGAACTTCC\|ACAGAGGGCAAAAAA\|ACTTCTGAGGCTCCA\|GTTATGCAAGGTCCC\|TGTGCTGTGATTGGCG\|TCAGCGTCCGACCCAG\|GATGCGCTCACCTCCC\|TCTTGGCCATGGTGCC\|GCGCCTGCGCACTGAG\|GAGACGGAGGAAGGTC\|GACTCCGGCCGGCTCT\|CCGCGTGCGTCCCGGGTTAC\|CCATCTGTCCTAATTCGGAT' > matched.fastq
grep太快了，R望尘莫及！
二代reads最短都有50bp，所以大家常用的比对工具都是不支持50bp以下的reads的比对的。
但是，在实际中，我们确实又有比对super short reads的需求。
So，我找到了如下方法来比对super short reads：
1.msa.sh，bbmap里的，好用，但是太慢。
2.bbduk.sh，bbmap里的，不支持输出sam文件，只能输出map和unmap的reads。
3.bbmap.sh，好用，超快，但是不支持mismatch，可惜。
4.seal.sh，没试过。
5.bowtie1，亲测可用，输出Sam文件。
