癌症样本全转录组数据的融合基因鉴定
前几期转录组周更学习分享了lncRNA和mRNA联合分析的一般套路和鉴定新lncRNA的基本流程，接下来的两周我会带大家一起学习之前一位老师对癌症样本全转录组数据进行融合基因和变异鉴定的推文 老程的全转录组，解决遇到的各种问题
本次教程使用数据为GSE145894，是一套胃癌的全转录组测序，测序平台为Illumina，方法为PE150。 至于mRNA-seq，如果能够处理全转录组数据，那么更常见的polyA富集建库的RNA-seq自然不在话下。
在这项研究中，我们建立了一种连续异种移植模型，通过在小鼠体内用5-氟尿嘧啶（5-Fu）长期治疗EBVaGC细胞株SNU719，来富集EBVaGC癌干细胞（CSCs）。
简单地说，给皮下注射SNU719的NOD/SCID小鼠每周用5-Fu治疗，直到异种移植物直径达到1.5厘米，然后将小鼠安乐死，用胶原消化法获得单细胞悬浮液。然后按上述方法将纯化的细胞传给经 5-Fu 处理的 NOD/SCID 小鼠。连续四代后，我们从经5-Fu处理的第四代异种移植中获得了SNU-4th细胞。
与亲代SNU719细胞相比，SNU-4th细胞具有干性特征。越来越多的证据证明，环状RNA（circRNA）在维持肿瘤表型方面发挥着重要作用。因此，我们通过RNA测序比较了用5-Fu和PBS处理的第四期异种移植物的circRNA表达，希望找出哪些circRNA参与了EBVaGC的干性调控。我们通过用5-氟尿嘧啶长期处理EBVaGC细胞株SNU719，建立了高度恶性的EBVaGC细胞株SNU-4th，然后用RNA测序法比较了用5-氟尿嘧啶和PBS处理的第四期异种移植物的circRNAs表达。
下载数据
kingfisher get -p PRJNA608689 -m ena-ascp ena-ftp prefetch aws-http
fastp
ls SRR111783* | cut -d "_" -f1 | uniq > ../srr_acc_list.txt
cat srr_acc_list.txt | while read id
do
fastp -i 0.raw/${id}_1.fastq.gz \
-o 1.qc/${id}_1.fastp.fq.gz \
-I 0.raw/${id}_2.fastq.gz \
-O 1.qc/${id}_2.fastp.fq.gz \
-l 36 -q 20 --compression=6 \
-R ${id} -h ${id}.html
done
我还在中山实习的时候曾老师就吐槽过我作为一个练习时长半年的生信菜鸟是时候学会使用自动化工具管理生信流程啦，所以后面的转录组专辑我可能会给大家带来如何使用snakemake管理流程，因为马上开学了没啥时间，这两周还是一个个脚本命令执行吧，大家再忍一忍orz
并且我自认为我的笔记一个个脚本整理成推文流程还是很完整的，虽然截图有点多，但起码如果大家真的跟我走下来工作路径、文件名称位置啥的不会成为障碍
结果文件：
运行log：
挑第一个html报表看看：
RNA-seq数据分析完全指北-02：fastq文件质量控制
- GC双峰
可以看到文件名、reads总数、读长是PE150，GC含量53%等信息。正常人类基因组的GC含量应该是43%左右，但是由于建库过程中不可避免的扩增偏倚等原因，最终的GC含量大概是50%左右。
这一点和我们fastp结果一致
这里有点问题，虽然总体GC含量基本正常，但是在理论分布的右侧还有一个GC峰，可能是rRNA，因为这次的示例数据用的是去除rRNA建库的，而这种建库方式一般很难完全将rRNA除去
可以看到基本所有的数据都有明显的GC双峰
fastqc -t 16 -o ./1.qc_html/ ./0.raw/SRR11178348_1.fastq.gz
fastp报表没有这个“per sequence GC content”，也看看我们的fastqc，一致
- 另外我们的fastqc报表和fastp报表都可以看出存在许多重复序列：
虽然对于RNA-seq来说，重复序列不可避免，但是高达80%的重复片段肯定也是有问题的，需要我们进一步对数据进行校正。（原推文单文件fastqc，我们这里也一样，sequence duplication level为1的只有百分之10多）
去除奇怪的RNA
1、GC双峰和重复序列来自哪里？
2、下载rRNA序列
jimmy曾在一篇推文中提到，去除rRNA可以去除GC双峰的右峰
可以发现rRNA多了许多
下载：
3、Hisat2构建索引并输出未比对的fastq序列
构建索引
hisat2-build -p 4 rRNA.fasta rRNA
输出没有比对到rRNA的序列
for i in {48..53}
do
hisat2 -x 2.filter_otherRNA/rRNA -1 0.raw/SRR111783${i}_1.fastq.gz -2 0.raw/SRR111783${i}_2.fastq.gz --un-conc-gz ./2.filter_otherRNA/SRR111783${i}_rmr_%.fq.gz -p 16 -S ./2.filter_otherRNA/SRR111783${i}.sam
done
-x 这个参数指定了Hisat2需要使用的参考基因组索引文件,这里是我们前面构建好的已知rRNA的索引
-1 -2 分别指定了配对的第一端和第二端测序数据的文件路径，${i} 被循环中的当前整数取代，用来构建不同的文件名
--un-conc-gz 指定了未比对成功的测序数据的输出路径，我们这里需要输出没有比对到rRNA的序列
%.fq.gz 部分会被Hisat2替换为输出文件的前缀，以便生成两个未比对的文件
-p 16: 这个参数指定了线程数，将在16个线程上并行运行Hisat2
-S 这个参数指定了比对结果的输出文件路径，生成SAM格式的比对结果
可以发现，匹配上rRNA的序列大部分都比对上超过一次
看看fastq文件size比较：
可以看到，去除rRNA序列之后，fastq文件大小都减少了大约百分之20，也可以通过查看log查看细节
再次质控看看，成功将右峰去除
fastqc -t 16 -o ./2.qc_html/ 2.filter_otherRNA/SRR11178348_rmr_1.fq.gz
比较去除rRNA之前的结果：
成功将右峰去除
ribodetector
曾老师给我推荐了一个自动去除rRNA的软件，不用自己去下载rRNA序列那么麻烦，我们这里也顺带使用比较一下
nohup ribodetector_cpu -t 10 -l 100 -i ./0.raw/SRR11178348_1.fastq.gz ./0.raw/SRR11178348_2.fastq.gz -e rrna --chunk_size 256 -r 2.ribodetector_test/reads.rrna.{1,2}.fq.gz -o 2.ribodetector_test/reads.nonrrna.{1,2}.fq.gz &> ribodetector_log &
查看去除rRNA后结果：
为什么没匹配上rRNA的序列文件比原文件大？
- 解压看看和原fq文件的序列数量关系：
原fq文件
$grep -c "^@" SRR11178348_rmr_1.fq
49273485
$grep -c "^@" SRR11178348_rmr_2.fq
49273485
ribodetector过滤掉rRNA的fq文件
$grep -c "^@" reads.nonrrna.1.fq
48605115
$grep -c "^@" reads.nonrrna.2.fq
48605115
发现序列还是比原文件少的，大小可能受压缩程度的影响
- 比较于使用hisat2比对从NCBI下载下来的rRNA fa文件去除结果：
使用hisat2比对从NCBI下载下来的rRNA fa文件去除后的fq文件
$grep -c "^@" SRR11178348_rmr_1.fq
42234765
$grep -c "^@" SRR11178348_rmr_2.fq
42234765
看起来使用hisat2比对从NCBI下载下来的rRNA fa文件去除了更多的序列（42234765< 48605115）
关于rebodetector，我在检索学习的时候发现TBtools的CJ老师也推过，这里就简单地介绍一下：
NAR：测序数据鉴别和去除rRNA序列利器RiboDetector
RiboDetector是一款用于从宏基因组、宏转录组、ncRNA和核糖体测序数据中准确而快速地检测和去除rRNA序列的软件（https://github.com/hzi-bifo/RiboDetector）。它是基于深度学习的BiLSTM（一款双向的循环神经网络架构）开发的软件。跟基于比对和隐马尔科夫模型的方法相比RiboDetector能抓取更长距离的序列特征，从而具有更好的准确性。文章中对比了目前常用的其他5款工具。在测试数据上，RiboDetector比其他软件的错误预测率低6到2000倍。另外它的CPU模式运行比目前最常用的软件有10倍左右的速度提升，而在GPU模式上实现了50倍左右的运行速度提升（见下图B-C）。文章中分析表明RiboDetector有很好的泛化能力，能预测发现新的rRNA（和数据库中已知rRNA相似性低于90%）序列。最后，在测试数据上它的假阳性预测序列没有显著的对某些功能的偏向性（没有GO term显著富集在假阳性预测的序列中）。
https://academic.oup.com/nar/article/50/10/e60/6533611
The presence of abundant rRNA can introduce substantial bias to transcriptome data. A number of protein coding genes contain some rRNA-like sequence segments (17–20), so their expression levels can be strongly overestimated unless the rRNA reads are removed. Furthermore, rRNA removal can greatly reduce the data size for downstream analysis and accelerate the entire workflow.
To facilitate the removal of rRNA sequences from large-scale sequencing data, a number of methods have been developed, including Meta-RNA (21), rRNASelector (22) and rRNAfilter in the RNA-QC-chain pipeline (23) (named RQC_rRNAfilter below), which use hidden Markov models (HMMs) trained on curated rRNA sequences; RiboPicker (24) and SortMeRNA (25), which are based on sequence alignment; and rRNAFilter (26), which applies an expectation-maximization algorithm to discriminate rRNA reads and non-rRNA reads based on their k-mer profiles. Some short-read aligners such as BWA (27) can also be used to remove rRNA reads.
我们这里使用的hisat2应该也属于： Some short-read aligners such as BWA (27) can also be used to remove rRNA reads. 测序进入大样本时代，长读长测序如何成为主导？
Most of these methods are either based on alignment algorithms that search for sequence homologies or based on a first-order HMM, which makes probabilistic predictions based on the assumption that the current state only depends on the state one step before.
Moreover, RNAseq datasets are usually large: one sample can comprise 10G bases and one dataset may have 100 such samples, and the existing methods take hours to process one such sample. Thus, those methods may require weeks to remove rRNA from an entire dataset.
- RiboDetector的目的和原理 ：RiboDetector是一种基于深度学习的方法，利用双向长短期记忆（BiLSTM）模型从长距离上下文中捕获模式，快速准确地检测rRNA序列。
- RiboDetector的训练和验证数据 ：RiboDetector使用来自SILVA数据库的rRNA序列和来自OMA数据库的编码序列（CDSs）作为训练和验证数据，通过聚类和过滤去除可能的rRNA污染。
- RiboDetector的性能评估 ：RiboDetector与其他基于比对或隐马尔可夫模型（HMMs）的方法在多个基准数据集上进行了比较，包括SILVA rRNA序列、16S rRNA基因扩增序列、OMA CDSs、病毒序列和ncRNA序列。结果显示，RiboDetector具有最低的误分类率、最快的运行时间、最小的输出大小和最低的偏差。
- RiboDetector的泛化能力 ：RiboDetector还展示了对与训练数据不同的rRNA序列的显著泛化能力，即使是与训练数据相差20%以上的rRNA序列，也能被有效地识别。而其他方法在这方面表现不佳，随着序列差异的增加，误分类率也随之增加。
当然去除了rRNA对全转录组的数据来说还不算完，去除rRNA后右峰消失基本上可以认为右峰是rRNA造成的，但是可以发现还存在一个小小的左峰
4、左峰是什么？
既然右峰是rRNA，那么左峰有没有可能是tRNA呢？
构建索引
hisat2-build -p 4 tRNA.fasta tRNA
输出没有比对到tRNA的序列
for i in {48..53}
do
hisat2 -x 2.filter_otherRNA/tRNA -1 2.filter_otherRNA/SRR111783${i}_rmr_1.fq.gz -2 2.filter_otherRNA/SRR111783${i}_rmr_2.fq.gz --un-conc-gz ./2.filter_otherRNA/SRR111783${i}_rmr_rmt_%.fq.gz -p 16 -S ./2.filter_otherRNA/SRR111783${i}_tRNA.sam
done
去除tRNA前后fastq文件大小比对，基本一样：
查看log文件：
和原推文一样，我们没有比对上tRNA
5.其他序列
其实还有一些结构性RNA需要去除，包括scRNA、SRP RNA还有Ribonuclease P RNA Component H1等，获得这些序列的方法类似，但是过程要更加繁琐一些，这里就不具体介绍了。我把整理好的otherRNA.fa（包括rRNA、tRNA和otherRNA）上传到了百度网盘，需要的读者可以自取。
上面的百度网盘连接已经失效了，我直接找曾老师要了otherRNA.fa，因为我个人也不清楚如何从NCBI下载囊括所有的otherRNA，希望有知道的老师同学可以在评论区告诉我【Q1】
这个otherRNA.fa文件我也会放在公众号后台，回复关键词【otherRNA】获取
查看fa文件：
和原推文rRNA数量差不多
那就接着往下做
for i in {48..53}
do
hisat2 -x 2.filter_otherRNA/other_RNA/otherRNA -1 2.filter_otherRNA/SRR111783${i}_rmr_1.fq.gz -2 2.filter_otherRNA/SRR111783${i}_rmr_2.fq.gz --un-conc-gz ./2.filter_otherRNA/other_RNA/SRR111783${i}_rmr_rmOther_%.fq.gz -p 16 -S ./2.filter_otherRNA/other_RNA/SRR111783${i}_tRNA.sam
done
fastq文件大小比对：
这里就不看具体的序列数量了，解压怪麻烦的
可以发现不同样本去除程度这次并不一样
原推文中又做了fastqc、multiqc并提到：
可以看到前5对fastq文件现在质量已经可以勉强使用了，但是最后一个文件仍然有很大的问题。这么奇怪的GC含量，会不会是有其他物种污染呢？
并通过这篇推文进探索是否存在一些其他物种的污染：
RNA-seq数据分析完全指北-04：创建本地blast库分析物种组成
但最终没有比对上其它物种信息
所以我们在这里也不深究这个问题了，因为我们的主要目的是走通并学习这个流程
去接头以及过滤
我们前面fastp其实已经走了这一步，只不过针对的是raw fastq文件
仅从第一个样本文件可以看到这个时候基本上已经没有什么接头序列了：
但是原推文推文中显示，去除otherRNA后第三个样本的reads1有超过30%的reads存在接头：
所以这里我们还是使用trim-galore走一下质控，fastqc看看
关于trim-galore转录组专辑前面多次用到，需要注意这里使用时参数的变化
这里随便粘一个之前使用的脚本代码：
cat $1 |while read id
do
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o ./ $fq1 $fq2
done
原推文这里使用的脚本代码：
trim_galore --length 35 --paired --retain_unpaired --cores 16 -o ./ ${id}_rm_1.fq.gz ${id}_rm_2.fq.gz
参数解读：
--paired This option performs length trimming of quality/adapter/RRBS trimmed reads for paired-end files. To pass the validation test, both sequences of a sequence pair are required to have a certain minimum length which is governed by the option --length (see above). If only one read passes this length threshold the other read can be rescued (see option --retain_unpaired). Using this option lets you discard too short read pairs without disturbing the sequence-by-sequence order of FastQ files which is required by many aligners.
PE测序需要指定这个参数
-j/--cores INT Number of cores to be used for trimming [default: 1].
## 多线程运行
-q/--quality <INT> Trim low-quality ends from reads in addition to adapter removal. Default Phred score: 20.
## 碱基质量小于多少的reads会被去除，一般选择Q20，严格的话也可以选用Q30
--phred33 Instructs Cutadapt to use ASCII+33 quality scores as Phred scores (Sanger/Illumina 1.9+ encoding) for quality trimming. Default: ON.
## 默认fastq文件使用phred33
--stringency <INT> Overlap with adapter sequence required to trim a sequence. Defaults to a very stringent setting of 1.
## 在3‘末端至少有几个连续碱基被认为是接头序列之后会被去除，默认是1个
这几个质控过滤参数有默认值，不设置也无伤大雅
--length <INT> Discard reads that became shorter than length INT because of either
quality or adapter trimming. A value of '0' effectively disables
this behaviour. Default: 20 bp.
For paired-end files, both reads of a read-pair need to be longer than
<INT> bp to be printed out to validated paired-end files (see option --paired).
If only one read became too short there is the possibility of keeping such
unpaired single-end reads (see --retain_unpaired). Default pair-cutoff: 20 bp.
## 在去除接头之后，长度小于这个值的reads将被舍弃。对于PE测序，一个pair两个reads的长度都要大于这个值。当然也可以保留其中超过阈值的一条而舍弃另一条。
在去除接头之后，长度小于这个值的reads将被舍弃。对于PE测序，一个pair两个reads的长度都要大于这个值。当然也可以保留其中超过阈值的一条而舍弃另一条。
这时就需要看看原推文脚本代码这里增加的一个参数 --retain_unpaired ：
--retain_unpaired If only one of the two paired-end reads became too short, the longer
read will be written to either '.unpaired_1.fq' or '.unpaired_2.fq'
output files. The length cutoff for unpaired single-end reads is
governed by the parameters -r1/--length_1 and -r2/--length_2. Default: OFF.
## 保留pair中通过length阈值的一条，舍弃另一条
原推文中加上这个参数做出来的结果：
SE只有不到10M，相对于总量2G的测序数据不值一提，所以可以舍弃，后续分析会比较简便
所以我们这里就直接不添加这个参数，要求一个pair两个reads的长度都要大于指定的 --length 参数
这个软链接的骚操作之前没用过
cat srr_acc_list.txt | while read id; do (echo "trim_galore --length 35 --paired --cores 2 -o 3.trimg/ 3.trimg/${id}_rmr_rmOther_1.fq.gz 3.trimg/${id}_rmr_rmOther_2.fq.gz "); done > trim.sh
再次质控：
fastqc -t 16 -o ./3.trimg_qc ./3.trimg/*.fq.gz
multiqc ./3.trimg_qc/*zip -o ./
因为前面我们在3.trimg中放了软链接，所以这里我们相当于是同时对trim-galore前后的fq文件质控
查看multiqc报表和生成日志：
我们这里主要看两个部分
- trim-galore前后接头情况
- trim-galore前一个fq文件有多少reads有接头序列，和fastp的报告不同吗
在记录日志中
trim-galore日志中记录的每个fq文件Reads with adapters一般是30%多
这里用fastp对去除了otherRNA的第一个样本QC看看：
fastp -i SRR11178348_rmr_rmOther_1.fq.gz -o SRR11178348_rmr_rmOther_1.fastp.fq.gz -I SRR11178348_rmr_rmOther_2.fq.gz -O SRR11178348_rmr_rmOther_2.fastp.fq.gz -l 35 -q 20 --compression=6 -R SRR11178348_rmr_rmOther -h SRR11178348_rmr_rmOther.html
从结果文件大小上看trim-galore和fastp差不多
The input has little adapter percentage (~0.566311%), probably it's trimmed before.
那么这里fastp的adapter percent和trim-galore记录中的Reads with adapters有什么不同呢，希望有知道的老师同学在评论区告诉我【Q2】
前面因为是全转录组测序的数据，相关过滤、分析、报表我们看的比较仔细，现在我们正式进入寻找融合基因和gatk找突变的学习
比对到基因组并定量
STAR全称Spliced Transcripts Alignment to a Reference，是大名鼎鼎的ENCODE计划的御用RNA-seq比对软件。STAR使用底层的C++语言编写，可以多核运行，具有极快的比对速度。相比于另外两款常用有参转录组比对软件TopHat和Hisat2具有更高的唯一比对率。与GATK良好的兼容使RNA-seq寻找基因突变更加容易。此外，10X的单细胞转录组上游软件Cellranger也基于STAR构建。
所以我们学习单细胞可以看到cellranger也是一个软件走比对、定量
STAR的基本流程包括两步： 1.构建基因组索引：这一步用户需要提供基因组参考序列（FASTA文件）和注释文件（GTF文件）。仅需要构建一次就可以用于以后的所有对比了。 2.将reads比对到基因组上。
关于参考基因组和注释文件的选择：
以及参考： 如何选择参考基因组和注释文件(https://blog.csdn.net/flashan_shensanceng/article/details/115705200) 转录小白日志丨关于参考基因组的那些事儿(http://www.frasergen.com/cn/info_173.aspx?itemid=1085) 提供primary_assembly，选它！提供dna_sm，选它！提供dna_sm.primary_assembly，务必选它选它！ gff3/gtf：全部的注释信息（转录分析就选它！）
根据我个人的使用经历，dna.toplevel和dna_sm.toplevel都用过
文件名组成主要涉及两种组装形式和三种重复序列处理方式：
转录小白日志丨关于参考基因组的那些事儿(http://www.frasergen.com/cn/info_173.aspx?itemid=1085)
参考这些资料，应该是dna_sm.primary_assembly最合适
https://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz
使用来自GENCODE的注释文件 gencode.v44.annotation.gtf
参考猪狗的参考基因组构建索引构建STAR索引
注意这里建索引就用到了gtf注释文件，后面比对时反而没再用
nohup STAR --runMode genomeGenerate --genomeDir ./star/human --genomeFastaFiles ../Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa --sjdbGTFfile ../gencode.v44.annotation_RmChr.gtf --sjdbOverhang 149 --runThreadN 12 &
fastq文件比对到基因组：
cat srr_acc_list.txt | while read id; do echo -n "STAR --runThreadN 12 "; echo -n "--genomeDir ~/refdata/index/star/human "; echo -n "--outSAMtype BAM SortedByCoordinate --outReadsUnmapped Fastx "; echo -n "--quantMode GeneCounts --readFilesCommand zcat --twopassMode Basic "; echo -n "--outFilterType BySJout --outFilterMultimapNmax 20 "; echo -n "--outFilterMismatchNmax 999 --outFilterMismatchNoverReadLmax 0.04 "; echo -n "--alignSJoverhangMin 8 --alignSJDBoverhangMin 1 "; echo -n "--chimSegmentMin 20 --chimJunctionOverhangMin 20 --chimOutJunctionFormat 1 "; echo -n "--alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 "; echo -n "--chimSegmentReadGapMax 0 --alignSJstitchMismatchNmax 0 -1 0 0 "; echo "--readFilesIn ~/fusion_and_mut/3.trimg/${id}_rmr_rmOther_1_val_1.fq.gz ~/fusion_and_mut/3.trimg/${id}_rmr_rmOther_2_val_2.fq.gz --outFileNamePrefix ${id}"; done > star.sh
关于上述命令，原推文中做出了各部分的中文解释，还画了个草图：
- 比对定量部分
STAR --runThreadN 12 # 12线程
--genomeDir ~/reference/linux/STAR/STAR_GRCh38_genecode_v33/ref_genome.fa.star.idx/ # 参考基因组索引所在位置
--outSAMtype BAM SortedByCoordinate # 输出经过坐标排序的BAM文件
--outReadsUnmapped Fastx # 输出没能比对到基因组上的序列，格式与输入文件相同
--quantMode GeneCounts TranscriptomeSAM # 输出基因的Read Count文件以及转录本定量的SAM文件
--readFilesCommand zcat # 输入的fastq文件经过gzip压缩
--twopassMode Basic # STAR特有，两次对比模式
--readFilesIn ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz # 输入文件的名称
--outFileNamePrefix ${id} # 输出文件的前缀
## 以下参数设置来自ENCODE官方，有些解释很难翻译成中文，参见下图
--outFilterMultimapNmax 20 # 如果一个读段被多重比对超过20次，则认为这个读段不能被比对到基因组
--outFilterMismatchNmax 999 # 每对读段允许错配999个碱基（相当于不过滤）
--outFilterMismatchNoverReadLmax 0.04 # 每对读段允许出现读长*4%的碱基错配，即PE150允许2*150*0.04=12个碱基错配
--alignIntronMin 20 # 内含子最短是20个碱基
--alignIntronMax 1000000 # 内含子最长是1000000个碱基
--alignMatesGapMax 1000000 # 一对读段之间最长距离为1000000个碱基
- 可变剪切部分
## 以下参数设置来自ENCODE官方，有些解释很难翻译成中文，参见下图
--outFilterType BySJout # 对junction进行过滤以减少错误
--alignSJoverhangMin 8 # 未注释过的junction的最少的overhang是8个碱基
--alignSJDBoverhangMin 1 # 注释过的junction的最少的overhang是1个碱基
## 其他参数
--alignSJstitchMismatchNmax 0 -1 0 0 # 允许剪切点错配的个数（-1代表无限制）四个数字分别代表(1)非经典；(2)GT/AG或CT/AC；(3)GC/AG或CT/GC(4)AT/AC或GT/AT
- 融合基因部分
--chimSegmentMin 20 # 每对嵌合读段较短的一端至少有20个碱基，即PE150允许280+20结构的融合基因
--chimOutJunctionFormat 1 # 输出的Chimeric.out.junction文件可直接用于融合基因
--chimSegmentReadGapMax 0 # 嵌合读段之间不允许空位
--chimJunctionOverhangMin 20 # 嵌合的junction的最少的overhang是20个碱基，为了过滤非常短的外显子，即连续剪切事件
参数众多，仍不是很好理解
这里先联系其它资料挑重点学习：
参考基因组下载和构建STAR索引（人源）(https://www.jianshu.com/p/9bdad4a4f98f)
STAR \
--runMode genomeGenerate \ #让STAR执行基因组索引的生成工作
--runThreadN 10 \ #构建运行使用的线程数
--genomeDir . \ #构建好的参考基因组存放的位置，最好是单独建立的一个文件夹，这里是/references/Homo_sapiens/Ensembl/GRCh38
--genomeFastaFiles ./Homo_sapiens.GRCh38.fa \ #fasta文件（参考基因组序列文件）
--sjdbGTFfile ./Homo_sapiens.GRCh38.99.gtf \# gtf文件（基因注释文件）
--sjdbOverhang 149 #读段长度: 后续回帖读段的长度, 如果读长是PE 100， 则该值设为100-1=99
--shdbOverhang 我们常见二代测序短读长为150bp，所以这里常设设置为149
生信软件 | STAR（测序序列与参考序列比对）(https://cloud.tencent.com/developer/article/1849258)
STAR比对的基本原理：
在代码的参数中我们可以看到一些之前没见过的术语，如overhang、junction等，可以参考这脸文章学习：
一文搞懂基因融合（gene fusion）的定义、产生机制及鉴定方法(https://zhuanlan.zhihu.com/p/84928559#tocbar--jig2lh)
在STAR比对工具原理及使用介绍这篇推文中，作者进一步说明了什么是overhang，以及如何区分可变剪接和基因融合事件（当overhang同时匹配上原基因和其它基因外显子区域时）
结果解读：
- log文件
看到这里，是不是觉得STAR不仅仅是在做比对，甚至还帮我们做完了质控！ 有这么多质控结果，那么，我们能不能用把所有质控结果综合起来呢？ 我去multiqc上查了查，multiqc居然也支持STAR的质控报告！！这就方便多了！
于是我们这里用multiqc汇总Log.final.out看看：
可以发现，最后一个样本很明显在比对参考基因组时大部分没有比对上 Unmapped：too short
联系前面在去除otherRNA后质控发现的GC含量分布存在的问题
值得注意
- Aligned.sortedByCoord.out.bam文件
比对上并按照坐标排好序的bam文件
Unmapped.out.mate文件是未能比对到基因组上的片段重新生成的FASTQ文件
- SJ.out.tab
剪切文件记录的是Splice junctions相关信息，STAR 软件将intron的最后一个或者开始第一个碱基作为可变剪接的开始或者结束
$head SRR11178348SJ.out.tab
1 10102 10178 2 2 1 1 0 54
1 15039 15795 2 2 1 0 4 11
1 16311 16606 2 2 1 11 0 53
1 16766 16857 2 2 1 0 64 52
1 17056 17232 2 2 1 33 291 73
1 17369 17605 2 2 1 32 39 56
1 17743 17914 2 2 1 15 135 68
1 18062 18267 2 2 1 0 89 75
1 18062 24737 2 2 1 0 8 68
1 18062 195262 2 2 1 0 8 68
文件以 tab 作为分隔符，一共有9列：
- 染色体号
- intron的第一个碱基
- intron的最后一个碱基
- 正负链，其中0=undefined，1=+，2=-
- intron motif内含子剪切模式：0：非经典；1：GT/AG；2：CT/AC；3：GC/AG；4：CT/GC；5：AT/AC；6：GT/AT；
- 在可变剪接数据库中是否存在，0=unannotated，1=annotated
- 跨过该剪接位点的唯一比对reads数目
- 跨过该剪接位点的多处比对reads数目
- 该剪接位点的最长overhang（比对到的读段距离剪切点最远延伸了多少）
- Chimeric.out.junction
Chimeric.out.junction文件包含了可能的融合基因以及circRNA比对结果，可以进一步使用STARfusion或其他用于circRNA鉴定的软件进行鉴定
5 95763507 - 5 95783964 - 2 2 2 SRR11178348.621558
95763508 68S81M-79p111M4425N38M 95783896 68M81S
7 86979161 + 7 86978469 + 0 0 1 SRR11178348.621607
86979127 1S34M112S 86978470 35S112M-62p149M
1 1257126 + 1 1257001 + -1 0 0 SRR11178348.621759 1257002 124M 1257002 122M
16 20821761 - 16 20828538 - 2 1 2 SRR11178348.621823
20821762 58S89M-82p134M3930N16M 20828480 58M89S
17 30465061 - 17 30465513 - 0 0 4 SRR11178348.622172
30465062 60S90M-91p149M 30465453 60M90S
11 57737607 - 11 57739261 - 2 0 2 SRR11178348.622230
57737608 27S61M167N62M 57738987 53M91N95M8p27M123S
17 5357158 - 17 5360281 - 0 0 4 SRR11178348.622256 5357159 101S48M2p150M 5360180 101M48S
11 77629900 + 11 77619605 + 2 0 2 SRR11178348.622277
77629848 1S52M96S 77619606 53S92M-11p9M2804N140M
4 153270337 - 4 153276131 - 2 0 0 SRR11178348.622282
153270338 130S20M-20p149M 153276001 130M20S
10 72715786 - 10 72715885 - -1 0 0 SRR11178348.622293
72715787 98M 72715788 97M
后面使用STAR-Fusion鉴定融合基因时会用到
一般来讲，此文件包含14列：
以第一行为例
- ReadsPerGene.out.tab
这个文件就是常见基因表达矩阵的来源，我们这里使用的是来自GENECODE的gtf注释文件
$head SRR11178348ReadsPerGene.out.tab
N_unmapped 3054008 3054008 3054008
N_multimapping 1476173 1476173 1476173
N_noFeature 10134300 32394512 10427757
N_ambiguous 2631155 36468 1292347
ENSG00000290825.1 0 0 0
ENSG00000223972.6 0 0 0
ENSG00000227232.5 47 0 47
ENSG00000278267.1 0 0 0
ENSG00000243485.5 0 0 0
ENSG00000284332.1 0 0 0
文件以 tab 作为分隔符，一共有4列：
1.基因ID
2.不区分正负链时的计数结果
3.比对到1st read strand的计数结果
4.比对到2st read strand的计数结果
所以，如果我们是没有链特异性的数据，直接提取第1-2列信息即可形成我们的表达矩阵
ReadsPerGene.out.tab文件包含了每个基因的count数，并且与htseq-count等方法得到的结果高度一致
相关基础转录组下游分析在RNA-seq数据分析完全指北-08：差异分析 有详细的代码和结果，从差异表达、富集分析搭配GSEA
转录组专辑前面已经介绍很多，感兴趣小伙伴也可直接参考言推文，我们不再赘述
此外，还存在STARgenome和STARpass1，这两个临时文件夹
之所以存在这两个临时文件夹，主要是我们使用了STAR的2-pass mode，在命令中的体现 --twopassMode Basic # STAR特有，两次对比模式 参数
可以发现更加灵敏的new junction
参考：
比对软件STAR的简单使用 (https://www.bioinfo-scrounger.com/archives/288/)
接下来就是应该拿我们这里产生的Chimeric.out.junction文件使用star-fusion进行融合基因的鉴定
但是我在进行的时候发现了一样的参数跑下来没有发现任何融合基因事件，当然这个问题目前我已经解决了，鉴于本文已经很长了，索性下次star-fusion和gatk找突变一块讲
以上就是本周转录组周更全部内容，留下了两个我不懂的问题，希望有老师同学帮我解答（Q1&Q2）
- 我个人也不清楚如何从NCBI下载囊括所有的otherRNA，比如有哪些结构RNA、如何检索
- fastp的adapter percent和trim-galore记录中的Reads with adapters有什么不同