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里的序列

?

1 2	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

?

1 2	grep --line-buffered AAAAAGGGGAGAGAGG whitelist.txt grep --line-number AAAAAGGGGAGAGAGG whitelist.txt

　　

批量运行

?

1 2 3 4 5 6 7

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

?

1	grep --ignore-case --line-number ttggaccaggatgggcaccacccgtttaagagc raw_gRNA | cut -f1 -d':' > empty.grep

　　

非常慢的python代码也贴出来吧：

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

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找到，最终却什么都没有。

?

1	zcat CRISPR_STD_CKDL230008134-1A_H3TGMDSX7_S1_L001_R1_001.fastq.gz | grep AAAGAACAGGTCATAA

于是就对CRISPR lib做自己的处理，UMI-tools提取barcode，然后自己做gRNA count矩阵。

最好对R1做一个trim：

?

1 2

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

 

?

1 2 3 4 5 6 7 8 9 10 11 12 13

# 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

 

参考：

• 软件使用 | UMI-tools产生单细胞表达矩阵

 

---

 

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的技术问题：

• How does Cell Ranger correct barcode sequencing errors? 
• What is a barcode whitelist? 

 

?

1 2 3

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

　

?

1 2 3 4 5 6 7

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，是否有很强的变异。

 

?

1 2 3 4

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代码

?

1

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文件。

 

参考：Good way to align very, very short reads?