背景
我們在進(jìn)行高通量測序后，會對篩選得到的差異基因進(jìn)行qpcr驗證。而在很多器官中，因為基因的特異性表達(dá)，有些外顯子存在組織特異性。這就會導(dǎo)致qpcr的結(jié)果可能與高通量得到的結(jié)論不同。為了將qpcr與高通量的結(jié)果更好地結(jié)合，我們可以考慮將高通量結(jié)果中qpcr靶向的序列進(jìn)行定量，這樣可以更好地與qpcr的結(jié)果交互驗證。
思路
1、使用blast確定primer比對到的序列，并獲得靶序列在基因組/轉(zhuǎn)錄組上的區(qū)間（參考鏈接、參考鏈接)）
2、使用samtools view根據(jù)上一步獲得的區(qū)間，將高通量數(shù)據(jù)中的靶序列的比對情況提取出來
3、使用featurecount對這部分bam文件進(jìn)行重新定量
方法
1. 安裝blast
#構(gòu)建conda環(huán)境
conda create -n blast
conda activate blast
#安裝
conda install blast
2. 構(gòu)建數(shù)據(jù)庫
#查看幫助文檔
blastn -help
makeblastdb -help
#構(gòu)建比對庫:小鼠的庫構(gòu)建得很快
nohup makeblastdb -in ~/my_genome_reference/mouse/mm10/Mus_musculus.GRCm38.dna.primary_assembly.fa -dbtype nucl -parse_seqids -out mouse_mm10 &
3. blast
test.fasta的內(nèi)容
>Trp53_primer1
CTCTCCCCCGCAAAAGAAAAA
>Trp53_primer2
CGGAACATCTCGAAGCGTTTA
>P21_primer1
CCTGGTGATGTCCGACCTG
>P21_primer2
CCATGAGCGCATCGCAATC
>MLKL_primer1
AATTGTACTCTGGGAAATTGCCA
>MLKL_primer2
TCTCCAAGATTCCGTCCACAG
>Cdkn2a_primer1
CGCAGGTTCTTGGTCACTGT
>cdkn2a_primer2
TGTTCACGAAAGCCAGAGCG
>telomerase_primer1
CGGTTTGTTTGGGTTTGGGTTTGGGTTTGGGTTTGGGTT
>telomerase_primer2
GGCTTGCCTTACCCTTACCCTTACCCTTACCCTTACCCT
>36B4_primer1
ACTGGTCTAGGACCCGAGAAG
>36B4_primer2
TCAATGGTGCCTCTGGAGATT
>Sod1_primer1
AACCAGTTGTGTTGTCAGGAC
>Sod1_primer2
CCACCATGTTTCTTAGAGTGAGG
>Sod2_primer1
CAGACCTGCCTTACGACTATGG
>Sod2_primer2
CTCGGTGGCGTTGAGATTGTT
比對：
#比對
blastn -task blastn -query test.fasta -out test.output -db mouse_mm10 -outfmt 7 -num_threads 4 -evalue 100
-max_target_seqs 2
#注意??！evalue是用來衡量比對顯著性的值，越小越顯著，但是我們用的是基因組做的庫，比對的primer可能是轉(zhuǎn)錄組，因此可能會有g(shù)ap，所以為了盡可能有hit，可以將它挑的大一點，我們最后根據(jù)匹配的堿基數(shù)目去鑒定特異性
#task 選擇blastn，默認(rèn)為megablast，只能找到相似度90%以上的序列！
#max_target_seqs 用于調(diào)節(jié)最多匹配的核酸序列數(shù)目
#當(dāng)然，最好使用RNA庫在做一遍
注意：這里得-outfmt可以用特定的格式展示輸出
-outfmt <String>
alignment view options:
0 = pairwise, #顯示配對情況
1 = query-anchored showing identities,
2 = query-anchored no identities,
3 = flat query-anchored, show identities,
4 = flat query-anchored, no identities,
5 = XML Blast output,
6 = tabular,
7 = tabular with comment lines, #帶有列名注釋
8 = Text ASN.1,
9 = Binary ASN.1
10 = Comma-separated values
4. 整理比對結(jié)果
library(stringr)
library(dplyr)
library(plyr)
x<-read.table("test.output")
#這里的SOD1和MLKL的匹配序列太離譜了，我們先刪除，之后手動檢查后做
x<-x%>%ddply(.,.(V1),function(x){x%>%arrange(V11)%>%.[1,]})
x$gene=str_split(x$V1,"_",simplify = T)[,1]
x<-x%>%filter(!gene%in%c("MLKL","Sod1"))
x<-x%>%ddply(.,.(gene),function(x){start=x$V9;end=x$V10;start=min(c(start,end));end=max(c(start,end));data.frame(chr=x$V2,start=start,end=end)})
x<-x[!duplicated(x$gene),]
x<-data.frame(position=paste(x$chr,":",x$start,"-",x$end,sep = ""))
write.table(x,"product_position.txt",row.names = F,col.names = F,quote = F)
#配合samtools view的參數(shù)格式，可以直接這樣打印出來
paste(x$position,collapse = " ")
#5:115561162-115561218 4:89294391-89294498 17:29098417-29098501 17:13008246-13008338 11:69589811-69590667
5. 將所有的bam文件提取子集并在該區(qū)間內(nèi)定量
-
使用samtools提取子集
mkdir qpcr_quantity&&cd qpcr_quantity #確定所有的bam文件的位置 find ../ -name *bam >bamfile.list #提取子集 samtools view -h bamfile 染色體：start-end cat bamfile.list |while read id;do samtools view -h $id 5:115561162-115561218 4:89294391-89294498 17:29098417-29098501 17:13008246-13008338 11:69589811-69590667 >${id}.sub;done mv `find ../ -name *sub ` ./ 使用featurecounts定量
featureCounts -T 10 -a ~/my_genome_reference/mouse/mm10/Mus_musculus.GRCm38.99.gtf -o read.count -p -B -C -t exon -g gene_name *bam.sub
6. 使用R對readcount做簡單統(tǒng)計
x<-read.table("../qpcr_quantity/read.count",header = T)
colnames(x)
all_rowsums<-x[,-c(1:6)]%>%rowSums()
x<-x[all_rowsums!=0,]
x<-x[,-c(2:6)]
x
因為做q的時候都要做內(nèi)參的，這里在標(biāo)準(zhǔn)化的時候，可以以內(nèi)參的count數(shù)為1，然后其它基因與內(nèi)參作比較（類似于q的標(biāo)準(zhǔn)化方法）
7. 小結(jié)
該方法的缺點：最佳匹配不一定是目標(biāo)序列
1、由于我使用的是基因組構(gòu)建的比對庫，其中不包含轉(zhuǎn)錄本的可變剪切信息，因此會導(dǎo)致某些引物比對質(zhì)量很差，如該案例中#補充，MLKL的reverse primer只有13bp能比對上，后來通過igv發(fā)現(xiàn)，其1-13在exon9上，其余的序列分布在別的exon上，這樣比對的結(jié)果是非常差的，因此后續(xù)要嘗試用轉(zhuǎn)錄組建庫
2、這種比對是把primer分別比對的，與primer blast還是有差別的，因此需要人工校驗比對的結(jié)果，比如primer1在2號染色體上，primer2可能在2號和4號染色體上均有最佳比對。
3、這里提供兩個下來refseq數(shù)據(jù)庫的方法：從RefSeq數(shù)據(jù)庫批量下載微生物基因組 - 知乎 (zhihu.com)、下載refseq序列_馬志遠(yuǎn)的生信筆記的博客-CSDN博客
4、登錄到NCBI的FTP的方式：如何下載NCBI refseq？ - 簡書 (jianshu.com)
#人的下載方式
#鼠的下載方式
分割符
基于轉(zhuǎn)錄本的定量方法
背景
如上所述，我們是基于基因組進(jìn)行qpcr靶序列的定量，但是這種方法有時候并不準(zhǔn)確，因此我們需要從轉(zhuǎn)錄本水平進(jìn)行數(shù)據(jù)重比對，然后定量，方法與基于基因組的方法基本吻合
實現(xiàn)方法
1、下載轉(zhuǎn)錄本數(shù)據(jù)
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/10090/109/GCF_000001635.27_GRCm39/GCF_000001635.27_GRCm39_rna.fna.gz
2、構(gòu)建索引及比對
#注意這里使用的轉(zhuǎn)錄本數(shù)據(jù)，按理講應(yīng)該使用基于基因組的比對方法（因為hisat會考慮可變剪切，轉(zhuǎn)錄本水平不應(yīng)該考慮這個。）但是為了方便，還是用hisat2哈！
#構(gòu)建索引
hisat2-build -p 10 ../GCF_000001635.27_GRCm39_rna.fna refseq_grcm39
#比對
hisat2 -p $threads1 --dta -x $index_file -1 $fq1 -2 $fq2 | samtools sort -@ $threads2 -o qpcr_quantity/$id.sorted.bam -
#構(gòu)建index
samtools index -@ $threads $bam
3、對qpcr引物進(jìn)行blast，確定在轉(zhuǎn)錄本 上的位置（參考上文2，3，4步）
這里注意：在使用refseq庫的時候，blast會默認(rèn)把region的版本號去掉，如XM_030245923.2會變成XM_030245923,所以要與bam的header比較一下
4、從比對好的bam中把靶序列區(qū)域的比對情況提取出來（參考上文第5步（1））
5、基于第3步的結(jié)果，手動做一個靶區(qū)域的gtf，用于featurecounts（搜了很多，沒有辦法從fasta直接變成gtf）
#示例文件
# chr gene start end region
# 1 NM_007475 36B4 503 559 NM_007475.5
# 2 NM_007669 P21 114 198 NM_007669.5
# 3 NM_009877 Cdkn2a 103 210 NM_009877.2
# 4 NM_011434 Sod1 179 295 NM_011434.2
# 5 NM_013671 Sod2 293 385 NM_013671.3
# 6 XM_030243820 MLKL 1267 1443 XM_030243820.2
# 7 XM_030245923 Trp53 1466 1529 XM_030245923.2
my_gtf<-data.frame(seqid=x$region,
source=".",
feature="exon",
start=x$start,
end=x$end,
score=".",
strand=".",
phase=".",
attributes=paste("gene_name"," ",'"',x$gene,'"',sep = ""))
#注意，一定是改成tab分割，不然計數(shù)的時候會報錯
write.table(my_gtf,"../qpcr_quantity/my_gtf_qpcr_quantity.gtf",row.names = F,col.names = F,quote = F,sep = "\t")
6、featurecount定量
featureCounts -T 10 -a my_gtf_qpcr_quantity.gtf -o read.count -p -B -C -t exon -g gene_name -f -F GTF *bam.sub
7、查看定量的結(jié)果，并讀入R做統(tǒng)計即可！
https://mytuchuangwhq2.oss-cn-qingdao.aliyuncs.com/image-20220422113235245.png
