背景
我們?cè)谶M(jìn)行高通量測(cè)序后薪铜，會(huì)對(duì)篩選得到的差異基因進(jìn)行qpcr驗(yàn)證漆弄。而在很多器官中脚草，因?yàn)榛虻奶禺愋员磉_(dá)酷誓，有些外顯子存在組織特異性披坏。這就會(huì)導(dǎo)致qpcr的結(jié)果可能與高通量得到的結(jié)論不同。為了將qpcr與高通量的結(jié)果更好地結(jié)合盐数，我們可以考慮將高通量結(jié)果中qpcr靶向的序列進(jìn)行定量棒拂，這樣可以更好地與qpcr的結(jié)果交互驗(yàn)證。
思路
1、使用blast確定primer比對(duì)到的序列帚屉，并獲得靶序列在基因組/轉(zhuǎn)錄組上的區(qū)間（參考鏈接谜诫、參考鏈接)）
2、使用samtools view根據(jù)上一步獲得的區(qū)間攻旦，將高通量數(shù)據(jù)中的靶序列的比對(duì)情況提取出來(lái)
3喻旷、使用featurecount對(duì)這部分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ù)庫(kù)
#查看幫助文檔
blastn -help
makeblastdb -help
#構(gòu)建比對(duì)庫(kù):小鼠的庫(kù)構(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
比對(duì)：
#比對(duì)
blastn -task blastn -query test.fasta -out test.output -db mouse_mm10 -outfmt 7 -num_threads 4 -evalue 100
-max_target_seqs 2
#注意！牢屋！evalue是用來(lái)衡量比對(duì)顯著性的值掰邢，越小越顯著，但是我們用的是基因組做的庫(kù)伟阔，比對(duì)的primer可能是轉(zhuǎn)錄組辣之，因此可能會(huì)有g(shù)ap，所以為了盡可能有hit皱炉，可以將它挑的大一點(diǎn)怀估，我們最后根據(jù)匹配的堿基數(shù)目去鑒定特異性
#task 選擇blastn，默認(rèn)為megablast合搅，只能找到相似度90%以上的序列多搀！
#max_target_seqs 用于調(diào)節(jié)最多匹配的核酸序列數(shù)目
#當(dāng)然，最好使用RNA庫(kù)在做一遍
注意：這里得-outfmt可以用特定的格式展示輸出
-outfmt <String>
alignment view options:
0 = pairwise, #顯示配對(duì)情況
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. 整理比對(duì)結(jié)果
library(stringr)
library(dplyr)
library(plyr)
x<-read.table("test.output")
#這里的SOD1和MLKL的匹配序列太離譜了灾部，我們先刪除康铭，之后手動(dòng)檢查后做
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ù)格式，可以直接這樣打印出來(lái)
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對(duì)readcount做簡(jiǎn)單統(tǒng)計(jì)
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
因?yàn)樽鰍的時(shí)候都要做內(nèi)參的赌髓，這里在標(biāo)準(zhǔn)化的時(shí)候从藤，可以以?xún)?nèi)參的count數(shù)為1，然后其它基因與內(nèi)參作比較（類(lèi)似于q的標(biāo)準(zhǔn)化方法）
7. 小結(jié)
該方法的缺點(diǎn)：最佳匹配不一定是目標(biāo)序列
1锁蠕、由于我使用的是基因組構(gòu)建的比對(duì)庫(kù)夷野，其中不包含轉(zhuǎn)錄本的可變剪切信息荣倾，因此會(huì)導(dǎo)致某些引物比對(duì)質(zhì)量很差，如該案例中#補(bǔ)充舌仍，MLKL的reverse primer只有13bp能比對(duì)上妒貌，后來(lái)通過(guò)igv發(fā)現(xiàn)灌曙，其1-13在exon9上，其余的序列分布在別的exon上推姻，這樣比對(duì)的結(jié)果是非常差的，因此后續(xù)要嘗試用轉(zhuǎn)錄組建庫(kù)
2藏古、這種比對(duì)是把primer分別比對(duì)的，與primer blast還是有差別的拧晕，因此需要人工校驗(yàn)比對(duì)的結(jié)果，比如primer1在2號(hào)染色體上厂捞，primer2可能在2號(hào)和4號(hào)染色體上均有最佳比對(duì)输玷。
3、這里提供兩個(gè)下來(lái)refseq數(shù)據(jù)庫(kù)的方法：從RefSeq數(shù)據(jù)庫(kù)批量下載微生物基因組 - 知乎 (zhihu.com)靡馁、下載refseq序列_馬志遠(yuǎn)的生信筆記的博客-CSDN博客
4欲鹏、登錄到NCBI的FTP的方式：如何下載NCBI refseq？ - 簡(jiǎn)書(shū) (jianshu.com)
#人的下載方式
#鼠的下載方式
分割符
基于轉(zhuǎn)錄本的定量方法
背景
如上所述臭墨，我們是基于基因組進(jìn)行qpcr靶序列的定量赔嚎，但是這種方法有時(shí)候并不準(zhǔn)確，因此我們需要從轉(zhuǎn)錄本水平進(jìn)行數(shù)據(jù)重比對(duì)胧弛，然后定量尤误，方法與基于基因組的方法基本吻合
實(shí)現(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)建索引及比對(duì)
#注意這里使用的轉(zhuǎn)錄本數(shù)據(jù)损晤，按理講應(yīng)該使用基于基因組的比對(duì)方法（因?yàn)閔isat會(huì)考慮可變剪切，轉(zhuǎn)錄本水平不應(yīng)該考慮這個(gè)红竭。）但是為了方便尤勋，還是用hisat2哈！
#構(gòu)建索引
hisat2-build -p 10 ../GCF_000001635.27_GRCm39_rna.fna refseq_grcm39
#比對(duì)
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茵宪、對(duì)qpcr引物進(jìn)行blast斥黑，確定在轉(zhuǎn)錄本 上的位置（參考上文2，3眉厨，4步）
這里注意：在使用refseq庫(kù)的時(shí)候锌奴，blast會(huì)默認(rèn)把region的版本號(hào)去掉，如XM_030245923.2會(huì)變成XM_030245923,所以要與bam的header比較一下
4憾股、從比對(duì)好的bam中把靶序列區(qū)域的比對(duì)情況提取出來(lái)（參考上文第5步（1））
5鹿蜀、基于第3步的結(jié)果，手動(dòng)做一個(gè)靶區(qū)域的gtf服球，用于featurecounts（搜了很多茴恰，沒(méi)有辦法從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分割斩熊，不然計(jì)數(shù)的時(shí)候會(huì)報(bào)錯(cuò)
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)計(jì)即可分冈！
https://mytuchuangwhq2.oss-cn-qingdao.aliyuncs.com/image-20220422113235245.png
