不同转录组流程结果到底该如何比较
最近学徒刚刚完成了RNA-seq的训练，就随机抽取了一个公共数据库项目给他作为作业，是，研究者通过CRISPR/Cas9对 nasopharyngeal carcinoma (NPC) 细胞系的p53引入了一个杂合突变 R280T ，然后看mRNA profiles of CNE2 (control) and KO CNE2 cells 的差异。
可以看到是6个测序样本
作者提供了表达矩阵：GSE130398_fpkm.txt.gz
第一个流程是 hisat2+featureCounts这个看我的视频就好：
第二个流程是salmon
第三个流程在文章里面，但是找不到文章，除非去询问作者，Xiangya Hospital, Central South University 的 zhenqi qin
首先查看突变是否引入成功
首先查看bam文件的头：
@PG ID:hisat2 PN:hisat2 VN:2.1.0 CL:"/project/home/lyang/miniconda2/bin/hisat2-align-s --wrapper basic-0 -p 10 -x /project/home/lyang/refdata/hisat/human/hg38/genome -S SRR8980083.sam -1 /tmp/3190.inpipe1 -2 /tmp/3190.inpipe2"
流程使用的是hg38参考基因组 , 简单搜索目标基因：https://www.ncbi.nlm.nih.gov/gene/7157
GRCh38.p13 (GCF_000001405.39) 17 NC_000017.11 (7668402..7687550, complement)
也可以在gtf文件搜索：
chr17 HAVANA gene 7661779 7687550 . - . gene_id "ENSG00000141510.17"; gene_type "protein_coding"; gene_name "TP53"; level 2; hgnc_id "HGNC:11998"; tag "overlapping_locus"; havana_gene "OTTHUMG00000162125.10";
批量提取 tp53 坐标的小bam：
ls *sort.bam | while read id;do (samtools view -b $id chr17:7661779-7687550 > ${id%%.*}_tp53.bam );done
ls *_tp53.bam |xargs -i samtools index {}
这个 杂合突变 R280T ，需要找到坐标，首先需要知道转录本，有一个网页工具可以继续转换，但是我忘记了。
假装作者是对的，他们的实验的确是引入了这个突变吧。本来都想发出去了，但是学徒凭运气找到了这个位置，给大家过目：
然后看相关系数
三种文件都准备好了：
首先看 salmon这样的无需比对的流程结果和 hisat2+featureCounts的差异
可以看到，同一处理组的样本在不同流程下面得到的表达量直接的相关性，是高于不同组的，符合逻辑！
但是单独查看同一个样本的不同流程的表达量，如下所示：
可以看到，还是有不少基因在不同流程表现差异非常显眼！那同样的，我们需要检查这些基因，简单看看5个差异最大的基因吧。
在salmon里面，可以看到第一个基因有3个转录本：
ENST00000445593.6_1 ENSG00000104341.16_2
ENST00000521545.6_1 ENSG00000104341.16_2
ENST00000517924.5_2 ENSG00000104341.16_2
同样的，salmon的这个样本的结果如下：
Name Length EffectiveLength TPM NumReads
ENST00000445593.6 3173 2867.291 0.114792 4.192
ENST00000521545.6 915 609.515 0.000000 0.000
ENST00000517924.5 618 313.557 0.000000 0.000
但是在hisat2+featureCounts的，6个样本都有表达量。
ENSG00000104341.16 chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8;chr8 97775057;97775057;97775776;97775849;97805353;97805353;97805353;97805353;97815304;97815328;97815328;97815328;97816058;97816058;97816058;97816058;97819140;97819140;97819140;97819140;97825058;97825058;97825058;97851397;97851397;97851397 97776108;97776108;97776108;97776108;97805464;97805464;97805464;97805464;97815401;97815401;97815401;97815401;97816180;97816180;97816180;97816180;97819164;97819238;97819238;97819238;97825153;97825153;97825153;97851474;97853013;97852598 +;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+;+ 3197 5453 5150 5346 4250 4775 5271
现在问题是如何判断这个基因是否有表达量，还是需要知道坐标：
chr8 HAVANA gene 97775057 97853013
使用命令去检查bam文件
samtools view SRR8980083.sort.bam chr8:97775057-97853013|wc
发现这个基因的区域的确是有比对成功的reads，这就是我们所说的表达量。
第一个结论，是salmon会漏掉基因的真实表达量。
那么salmon是否会凭空捏造基因的表达量，所以我们还需要检查
chr3 HAVANA transcript 96617188 96618236 . - . gene_id "ENSG00000269028.3"; transcript_id "ENST00000600213.3"; gene_type "protein_coding"; gene_name "MTRNR2L12"; transcript_type "protein_coding"; transcript_name "MTRNR2L12-201"; level 2; protein_id "ENSP00000468991.1"; transcript_support_level "NA"; hgnc_id "HGNC:37169"; tag "not_best_in_genome_evidence"; tag "basic"; tag "appris_principal_1"; havana_gene "OTTHUMG00000175726.3"; havana_transcript "OTTHUMT00000430905.3";
使用命令去检查bam文件
samtools view SRR8980083.sort.bam chr3:96617188-96618236 |wc
的确是没有表达量的，那么为什么salmon会捏造这个表达量呢？