ixxmu
commented
9 months ago 鉴定lncRNA流程全套代码整理 by 生信菜鸟团
前两期周更我们通过一篇文章的复现整理了mRNA和lncRNA分析基本流程，但并没有涉及新lncRNA的鉴定，本周的推文本质上是我个人学习鉴定lncRNA的全套流程笔记，整合了我们公众号往期的资源，对代码进行了勘误更新，内容非常详实。
读前须知：本篇推文并非教学笔记，而是学习笔记，所以我并不会提前和大家说哪些地方有哪些坑或者在推文中排除掉走过的岔路，而是给大家看我每步的结果，带大家跳到坑里去再一块跳出来。
我们首先从[「生信技能树」LncRNA数据分析实战]这个公开视频入手，尝试跟着曾老师走完这个项目
「生信技能树」LncRNA数据分析实战(https://www.bilibili.com/video/BV1Zg4y187ff?p=8&vd_source=852ec8cbb4975dabedb5d1f798b80c2a)
学徒第5月-02-LncRNA实战(https://mubu.com/doc/ISk-Ev1tg)
下游分析
背景知识
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE59245
Long noncoding RNAs expression patterns associated with chemo response to cisplatin based chemotherapy in lung squamous cell carcinoma patients
Non-coding RNA profiling by array
肺鳞状细胞癌（SCC）患者对顺铂化疗的反应存在很大的变异性。LncRNAs是一种潜在的新型预测标记物，可以识别从化疗中受益的患者亚组，对治疗指导具有重要价值。
使用对接受顺铂化疗的晚期肺SCC患者的部分反应（PR）肿瘤与进行性疾病（PD）肿瘤的微阵列分析来鉴定差异表达的lncRNA，并通过定量实时PCR（qPCR）进行验证。结果：与PD样本相比，953个lncRNA持续过度调节，749个lncRNAs
曾老师在视频中首先介绍的是下游的分析，因为这个大家通过常规转录组的学习已经比较熟悉了
这部分我们就不老生常谈，需要补全这部分知识的同学可以参考我们前面的推文初探mRNA、lncRNA联合分析之下游，或者自己去看看这个视频前面一部分，重难点主要是id注释，这一点可以通过曾老师的AnnoProbe包来帮助进行，这篇推文初探mRNA、lncRNA联合分析之下游中我主要用的是biomaRt包
我们重点还是下面这部分
上游pipeline鉴定lncRNA
背景知识
这部分视频也有相关推文与之对应：
主要是复现这个项目：
Transcriptome Analysis Suggests the Roles of Long Intergenic Non-coding RNAs in the Growth Performance of Weaned Piglets
就是重新下载一个公共数据，然后进行新lncRNA鉴定和注释，分析部分主要是分成两个大块，首先是hisat2+stringtie流程，然后是组装好的gtf文件的后，细致的进行新lncRNA鉴定和注释
测序数据获取
GSE65983
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE65983
https://www.ncbi.nlm.nih.gov/bioproject/PRJNA275632
https://sra-explorer.info/
这个项目GSE搜不到 PRJ可以搜到
去NCBI 下载accession list，这个后面整理文件感觉还是挺好用的
https://www.ncbi.nlm.nih.gov/Traces/study/?query_key=4&WebEnv=MCID_64c1e6edad6bc82565cdac33&o=acc_s%3Aa
下载参考基因组构建hisat2索引，参考：
fastp
上游的前面一部分也可以参考之前推文初探mRNA、lncRNA联合分析之上游
-i -o -I -O 分别对应一个样本两个fq文件的输入输出
看看帮助文档
-l, --length_required reads shorter than length_required will be discarded, default is 15. (int [=15])
--length_limit reads longer than length_limit will be discarded, default 0 means no limitation. (int [=0])
-q, --qualified_quality_phred the quality value that a base is qualified. Default 15 means phred quality >=Q15 is qualified. (int [=15])
-z, --compression compression level for gzip output (1 ~ 9). 1 is fastest, 9 is smallest, default is 4. (int [=4])
-R, --report_title should be quoted with ' or ", default is "fastp report" (string [=fastp report])
-h, --html the html format report file name (string [=fastp.html])
以第一个样本为例看看fastp报表：
参考：
[Fastp：过滤二代测序数据]https://www.jianshu.com/p/6da36135e2d1
hisat2比对
上游前面和bulk mRNA差不多
异常退出：
一直给我返回samtools Usage 感觉只剩软件版本问题？
之前做bulk mRNA的版本：
现在：
我记得我装的是新的，这是为什么？想起来了：
我在之前安装trinity时记录到的情况
所以这后面如果使用trinity需要特别注意
重装：
重装samtools后，单独测试第一个样本：
stringtie 定量
stringtie -p 4 -G $gtf \
-o 4.stringtie_gtfs/$sample.gtf \
-l $sample 3.hisat2_bams/$sample.bam
目前这一步用的还是Assemble RNA-Seq alignments into potential transcripts.
-p number of threads (CPUs) to use (default: 1)
-G reference annotation to use for guiding the assembly process (GTF/GFF3)
-o output path/file name for the assembled transcripts GTF (default: stdout)
-l name prefix for output transcripts (default: STRG)
到这里，我们就获得了stringtie对每个样本组装到的gtf文件
可以发现如果直接基于bam文件定量，给到的 -G 参数就是一开始的参考基因组gtf文件，而不是需要多走一步组装后merge的gtf文件
而我们这里就是要鉴定新的lncRNA
合并gtf后鉴定新的lncRNA
注意这一步stringtie merge 使用的-G参数 gtf文件仍为下载下来的已有参考，只有定量时使用的gtf文件是组装后merge的gtf文件
关于为什么需要merge后的组装gtf文件，stingtie相比featureCount能鉴定新转录本：
有关stringtie的使用和解读可以参考前面bulk mRNAseq这两篇
使用TPM/FPKM/RPKM进行差异分析真的可以消除系统误差吗？
stringtie （组装-合并-定量）
Stringtie和FeatureCount都是用于转录组定量的工具，但它们使用不同的定量策略。
在FeatureCount的情况下，它 直接将读数映射到预先存在的注释（如参考基因组或基因注释），并将读数分配给相应的特征 ，然后可以用于量化。
另一方面，Stringtie 使用从头开始的方法从读取中组装转录本并生成新的注释文件。然后使用该组装的转录组来估计来自同一组读数的表达水平。 进行合并步骤以合并多个样本并生成一致转录组，该转录组可用于定量所有样本中的读数，从而提高准确性和再现性。
这两种方法各有优缺点。FeatureCount的分析速度更快、更直接，但它只依赖于预先存在的注释，可能会错过以前没有注释的新的或替代的转录本。Stringtie速度较慢，但可以捕获新的或替代的转录物，并提供了对转录物丰度的改进评估，尤其是对于低表达或重叠的转录物。
当然stringtie其实也可以直接定量，如果不需要发现新的转录本和基因的话，可直接基于 bam 文件走定量步骤，我们这里还是使用从头开的方法学习一下这个软件。从头开始的方法和软件具体命令参数可参考上面给出的文章。
我们这里使用的为 组装-合并-gffcompare鉴定和筛选新的lncRNA
我们这里使用110和99看看：
可以看到，使用不同版本的gtf文件，对结果的影响很明显哦，鉴定新的lncRNA的话，也许2016年你认为是新的lncRNA，到2018年就已经被记录在册啦。
关于gffcompare的使用，参考：
比较不同样本的转录本定量信息需要先将转录本信息储存为相同的格式，一般组装软件的输出结果都是gtf或gff。由于在组装的过程中产生了大量的新的转录本信息，而我们仅通过肉眼观察其唯一的注释信息----染色体上的起始位置，很显然无法阐明其中蕴含的生物学意义，因此我们需要将它们与已知的转录本注释文件---annotation.gtf进行比较，将新得到的转录本与注释好的转录本之间建立联系，这样可以让我们更好地发现新的转录本。而gffcompare就是做的这个工作，由于它是基于cufflinks的一个附件cuffcompare开发的，因此很多原理及输出文件的格式也与cuffcompare类似。
我的思考：
为什么stingtie组装时使用的基因组参考注释gtf文件和gffcompare使用的基因组注释文件一样，还能找到新转录本？
StingTie是一个用于转录组组装的软件工具，它将测序数据映射到参考基因组上，并尝试重建转录本信息。在StingTie组装过程中，它会生成一个结果文件，通常是以gtf或gff格式保存的转录本注释信息。
为了发现新的转录本，StingTie将生成的结果文件与已知的转录本注释文件（例如annotation.gtf）进行比较，并将新得到的转录本与已知的转录本建立联系。这个比较的过程是通过使用gffcompare工具来实现的。
虽然StingTie使用的基因组参考注释文件和gffcompare使用的基因组注释文件可能是相同的文件，但在具体的操作中，它们扮演了不同的角色。
StringTie使用基因组参考注释文件来辅助转录组组装过程，提供基因和转录本的位置信息等，但它并不能完全包含所有已知的转录本信息，也不能涵盖特定样本的表达情况。
而gffcompare则是根据StingTie生成的结果文件和已知的转录本注释文件进行比较，从而找到新的转录本。它能够对比两个转录本集合，并确定新的转录本、已知转录本的变体以及已知转录本的完整性等信息。
因此，尽管两者使用的基因组注释文件可能一样，但StingTie和gffcompare从不同的角度对待这个注释文件，通过组装过程和比较过程的不同，能够找到新的转录本。
解读gffcompare结果文件：
输出文件六个,前四个文件可以指定保存位置，后两个文件是跟输入的gtf文件保存在一个位置，并且都是以-o提供的前缀开头的
- gffcmp.annotated.gtf：存储的是StringTie组装的转录本与注释文件内的转录本的差别信息，通过class_code来表示
如果只有1个新组装的gtf与参考gtf比较，则结果文件为gffcmp.annotated.gtf，包含新组装gtf文件里所有feature的注释结果；
如果有多个新组装的gtf与参考gtf比较，结果文件为gffcmp.combined.gtf，将所有新组装gtf的注释结果合并到一起
class code分类
= : 预测转录本与参考转录本拥有完全相同的内含子
c : 预测转录本包含在参考转录本中 /
j : 预测的转录本与参考转录本共享至少一个剪切位点，可能是潜在的新型isoform
i ： 预测的转录子完全落入参考内含子中
o ： 预测的转录本的外显子与参考转录本的外显子有重叠
u ： 与参考转录本相比，预测的转录本是在基因间区
x : 预测的转录本的外显子与参考转录本重合但是在相反的链上 /
e ： 预测的单外显子转录本与参考转录本至少重合10bp的参考内含子长度，有可能是pre-mRNA
p ： 预测的转录本的参考转录本附近2kb的距离内，可能是聚合酶滑动产生的片段
r ： 预测的转录本有50%以上的碱基与重复序列重合
s : 预测的转录本内含子与参考转录本重合但是在相反的链上(likely a mapping error)
class_code分类的具体含义: “=” 代码表示此预测转录本与注释基因的所有内含子完全吻合，但它们在第一外显子(first exon)的起始端或最后外显子(last exon)的末端可能有差别。然而，这并不影响将“=”类重建转录本判定为已注释转录本。又如，转录本标有 “j” 类别代码，表明此转录本至少有一个内含子与已注释基因的内含子相同，而其他位置可能不同，据此可推断此类转录本可能是注释基因的一个新异构体(novel isoform)。另外 “i,o,u,x” 的分类符合lncRNA的特征，可用于lncRNA的识别过程。因此，“i，j，o，u，x”这5类转录本表示可能是新的转录本，符合lncRNA的要求，保留作为后续分析。
- gffcmp.stats：文件存储比对结果的准确性和预测率
数据总结和准确性评估（包括Base、Exon、Intron、Intron chain、Transcript和Locus共6个水平的敏感性Sensitivity和准确性Precision评估）
gffcompare常用于比较从头组装的转录本与已知转录本，从而评估从头组装pipeline的性能，包括敏感性Sensitivity和准确性Precision两方面，从6个feature水平进行评估（Base level 、Exon level 、Intron level 、Intron chain level 、Transcript level 、Locus level）
敏感性Sensitivity=TP/(TP+FN) （对原真假而言）
准确性Precision=TP/(TP+FP) （对做出的判断而言，区别于特异性：正确判断原错误为错误，同敏感性，对原真假而言）
TP：true positives，真阳性，表示新组装的feature与参考的feature一致
FN：false negatives，假阴性，表示参考的feature未出现在新组装的feature中
FP：false positives，假阳性，表示新组装的feature在参考feature中未出现
FP+TP代表新组装feature的总数
loci基因组 transcript转录本
- gffcmp.loci
- gffcompare.tracking 该文件记录了所有样本中转录本的匹配情况，特别是对于多个新组装gtf的文件
- gffcmp..refmap
这个文件包含四列信息，第一列ref_gene_id是gene symbol ，无symbol的给出的是ensemble的gene id; 第二列ref_id是指ensemble的transcript id; 第三列class_code 是“=”和“c”;第四列是cuff_id_list。这个文件指组装后与参考基因组几乎完全匹配的转录本
对于每个新组装的gtf文件，都会产生一个refmap文件。内容为对于每个参考转录本，哪个query转录本完全或部分匹配到该转录本上，有4列，分别为：参考gtf中的基因名/基因ID；参考gtf中的转录本ID；匹配类型；详细匹配信息
- gffcpm..tmap 包含了转录本的定量信息，如cov，FPKM等，可用于定量或筛选新转录本(ref_gene_id,ref_id,class_code,qry_gene_id,qry_id,num_exons,FPKM,TPM,cov,len,major_iso_id,ref_match_len)
对于每个新组装的gtf文件，都会产生一个tmap文件，该文件可用于后续过滤转录本。内容为对于每个新组装gtf中的转录本，哪条参考转录本与其匹配度最高，一般有12列：参考基因名/基因ID；参考转录本ID；匹配类型；新组装基因ID；新组装转录本ID；新组装转录本外显子数；FPKM；TPM；Coverage；Length；新组装gtf中该基因的主剪切本；新组装转录本与参考转录本匹配的长度
从新组装的gtf文件中提取特定class code的转录本进行后续筛选
视频课对应的后续流程：
这篇文章想找intergenic （LINC）所有只选class code="u"
但是由于后面的视频课丢失，也无对应的推文：
所以从这里开始，我们参考另一个项目代码，继续我们的流程：
Gffcompare 获取转录本组装情况
这个我们前面根据视频课两个推文已经获得
关键文件：以release99版本为例
后面都是以最新版110为例
step1：保留指定class_code的transcripts
为跟进流程，我们这里就不再跟视频课对应项目要求，而是使用现在这个参考项目要求，比如取class_code: u x j i o 而不仅仅是 u
过滤，只保留class_code="u","x","i","j","o"的 transcripts
获取剩余transcripts的exons位置信息，提取序列并组装成转录本序列
获取剩余的transcripts的ID
整理transcripts得到对应gtf
把filter1_transcript.gtf中的class_code "=" 替换为L， 去除剩余为去除的class_code "="
前面只取了class_code: u x j i o 怎么会有 =
最终（43475-609350，如一个transcript包含多个exon片段）
step2: 根据长度进行过滤
过滤，只保留exon>1并且长度>200bp的transcripts
筛选tmap文件
感觉没去多少
transcripts的exon组成gtf，根据exon位置信息提取基因组序列，组装成转录本序列：
transcripts的exon组成gtf
根据exon位置信息提取基因组序列，组装成转录本序列
这里gffread需要另外装一下
还怪好嘞自动帮我保留了索引文件
各transcript的fa：
step3：转录本编码能力预测
转录本编码能力预测，主要是4个软件，需要取交集：
参考:
CPAT
https://cpat.readthedocs.io/en/latest/
User need to download prebuilt logit model and hexamer table (https://sourceforge.net/projects/rna-cpat/files/v1.2.2/prebuilt_model/) for human, mouse, zebrafish and fly. For other species, we provide scripts to build these models (see below).
https://cpat.readthedocs.io/en/latest/#build-hexamer-table
https://cpat.readthedocs.io/en/latest/#build-logit-model
Build hexamer table
Example:
make_hexamer_tab.py -c Human_coding_transcripts_CDS.fa -n Human_noncoding_transcripts_RNA.fa >Human_Hexamer.tsv
-c CODING_FILE, **--cod=**CODING_FILE
Coding sequence (must be CDS without UTR, i.e. from start coden to stop coden) in fasta format. User can get CDS sequence of a bed file using UCSC table browser(http://genome.ucsc.edu/cgi-bin/hgTables?org=Human&db=hg19&hgsid=289407045&hgta_doMainPage=1) -n NONCODING_FILE, **--noncod=**NONCODING_FILE Noncoding sequences in fasta format
-n NONCODING_FILE, **--noncod=**NONCODING_FILE
Noncoding sequences in fasta format
看看前面下载的fa文件：
如何判断存不存在UTR？
于是重新下载CDS fa文件
make_hexamer_tab.py -c ~/refdata/Sus_scrofa.Sscrofa11.1.cds.all.fa -n ../5.lncRNA/gffcomp110_filter2_transcript_exon.fa > Sus_scrofa_Hexamer.tsv
Build Logit model
Example:
make_logitModel.py -x Human_Hexamer.tsv -c Human_coding_transcripts_mRNA.fa -n Human_noncoding_transcripts_RNA.fa -o Human
-c CODING_FILE, **--cgene=**CODING_FILE
Genomic sequnences of protein-coding RNAs in FASTA (https://en.wikipedia.org/wiki/FASTA_format) or standard 12-column BED (https://genome.ucsc.edu/FAQ/FAQformat.html#format1) format.
CDS vs mRNA(https://www.differencebetween.com/difference-between-cds-and-cdna/):
下载mRNA/cDNA fa文件
make_logitModel.py -x Sus_scrofa_Hexamer.tsv -c ~/refdata/Sus_scrofa.Sscrofa11.1.cdna.all.fa -n ../5.lncRNA/gffcomp110_filter2_transcript_exon.fa -o Sus_scrofa
Run CPAT on local computer
cpat.py -x Sus_scrofa_Hexamer.tsv --antisense -d Sus_scrofa.logit.RData --top-orf=5 -g ../5.lncRNA/gffcomp110_filter2_transcript_exon.fa -o cpat_output
--antisense
Also search for ORFs from the anti-sense strand.Sense strand (or coding strand) is DNA strand that carries the translatable code in the 5′ to 3′ direction. default=False (i.e. only search for ORFs from the sense strand)
--top-orf=N_TOP_ORF
Number of ORF candidates. Many RNAs have dozens of possible ORFs, in most cases, the real ORF is ranked (by size) in the top several. To increase speed, we do not need to calculate “Fickett score”, “Hexamer score” and “coding probability” for all of them. default=5
CPC2
CPC2为CPC 的升级版，发布于2017 年，是目前最新的lncRNA 鉴定工具，也代表着lncRNA 鉴定的最新研究进展。
在经过大量的特征选择后，CPC2 最终的特征主要包括四条：最长ORF 长度，ORF 的完整性，Fickett 分数以及等电点 (isoelectric point, pI)[39,40]。其中等电点特征主要是通过将最长ORF 翻译为氨基酸序列，而后根据氨基酸等电点这一理化性质计算而得。与大多lncRNA 鉴定工具相同，CPC2 也使用了支持向量机来构建分类器。
# 创建python2环境conda create -n py2test python=2.7
安装biopython
conda install biopython=1.70
安装CPC2
wget https://github.com/biocoder/CPC2/archive/refs/heads/master.zip
unzip master.zip
cd ~/CPC2-master/libs/libsvm/libsvm-3.18
make clean && make
lncFinder
LncFinder是一种新的lncRNA识别工具。基于六聚体的对数距离，多尺度结构信息和从快速离散傅立叶变换获得的理化特征。
为了确定最佳分类器，使用10倍交叉验证对五种广泛使用的机器学习算法进行了验证：逻辑回归，支持向量机（SVM），随机森林，极限学习机器和深度学习。
最终选择SVM作为LncFinder的分类器。经过全面的功能选择和模型验证方案的评估，LncFinder在多个物种上的表现优于几种最先进的工具。用户可以轻松，高效地使用新的数据集或不同的机器学习算法对LncFinder进行重新训练。
https://cran.r-project.org/web/packages/LncFinder/index.html
Example：
# Prediction without secondary structure-based features
library(LncFinder)
demo_DNA.seq <- seqinr::read.fasta("~/lncRNA_project/07.identification/step2/filter2_transcript_exon.fa")
Seqs <- demo_DNA.seq
### The first parameter of lnc_finder() indicates the unevaluated
### sequences; parameter "SS.features = FALSE" means predict sequences
### without secondary structure-based features; parameter
### 'format = "DNA"' means the format of input sequences is DNA;
### parameters "frequencies.file" and "svm.model" indicate predicting
### sequences with which model; parameter "parallel.cores" means the
### number of cores to conduct parallel computing, and "-1" means
### using all available cores.
result_1 <- LncFinder::lnc_finder(Seqs,
SS.features = FALSE,
format = "DNA",
frequencies.file = "human",
svm.model = "human",
parallel.cores = 1)
write.table (result_1, file ="~/lncRNA_project/07.identification/step3/lncFinder/lncFinder_result.txt", sep ="\t",row.names =TRUE, col.names =TRUE,quote =FALSE)
lnc_finder（）的第一个参数表示未评估的序列；参数“SS.features=FALSE”表示预测没有基于二级结构特征的序列；参数'format=“DNA”'表示输入序列的格式为DNA；参数“frequencs.file”和“svm.model”表示用哪个模型预测序列；参数“parallel.cores”表示进行并行计算的核数，“-1”表示使用所有可用核。
我们这里是猪的需要自己构建这两个文件
需要传入lncRNA的fa文件，看起来lncFinder本质上是一个分类器，基于机器学习，需要有自己原有的训练数据集（已知lncRNA）？那这个文件去哪找呢？
mRNA需要：
https://anaconda.org/bioconda/viennarna
Returns data.frame. The first row is RNA sequence; the second row is Dot-Bracket Notation of secondary structure sequences; the last row is minimum free energy (MFE).
这一步耗时过长，遂跳过lncFinder的使用
PLEK
PLEK软件通过序列的kmer构成来区分编码和非编码转录本，不需要通过比对来完成，所以运行速度较快，同时其性能受到测序错误的影响的概率较低，比较稳定。
wget https://sourceforge.net/projects/plek/files/PLEK.1.2.tar.gz
tar xzvf PLEK.1.2.tar.gz
这个软件有点奇怪，难道只能放在工作目录下？
删除前面安装的目录和环境变量信息，重新在工作目录下安装，成功运行
看起来时脚本里遇到报错了被异常捕获到
但是还是有一个结果文件：
完整，那就拿着接着走吧
怀疑是-out参数要给一个文件名但我给的是目录
第一列代表该转录本为coding还是non-coding, 第二列为打分值，打分值大于0为coding, 小于零为non-coding, 第三列为fasta文件中的序列标识符
less _result | grep -w 'Non-coding' |awk '{print $3}' |sed 's/>//g' > plek_nc_id.txt
3个软件取交集，根据获取到的lncRNA预测的id，获取后续所需gtf和fasta文件
根据获取到的lncRNA预测的id，获取后续所需gtf和fasta文件（原推文无，此处为我自行查阅资料和参考其他步骤编写）：
提取对应的fasta序列文件：
（这一步其实可以融入到获取对应的gtf文件后面
参考前面step2：
transcripts的exon组成gtf，根据exon位置信息提取基因组序列，组装成转录本序列
这里只是因为一开始我觉得gtf文件没有必要了就没提，单独提取的fasta。其实gtf文件也需要一路根据过滤步骤过滤保存，在进行定量以便过滤掉低表达量的lncRNA是必要的）
提取对应的gtf文件：
此处注意文件命名，带有exon字样的gtf文件为提取后保留exon的gtf文件，是用来写作fa序列文件的
后面如果通过这个途径获取对应的fa文件就需要提取gtf文件中exon部分，通过gffread来写
step4：比对到Pfam据库
使用Transeq 将转录本序列翻译为6个可能的蛋白序列
nohup transeq filter3_by_noncoding_exon.fa filter4_protein.fa -frame=6 &
比对到[Pfam数据库]http://pfam.xfam.org/，过滤 (E-value < 1e-5)
Pfam数据库是一系列蛋白质家族的集合，其中每一个蛋白家族都以多序列比对和隐马尔科夫模型的形式来表示
http://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/
https://www.ebi.ac.uk/Tools/pfa/pfamscan/
参考：
使用pfam-scan进行预测(https://www.jianshu.com/p/fb3bd3de1c38)
Pfam数据库蛋白编码能力预测(http://www.biotrainee.com/thread-1212-1-1.html)
Pfam包括两个质量级别的家族数据库：Pfam-A和Pfam-B。Pfam-A来自基础序列数据库Pfamseq，是根据最新的UniProtKB数据建立的，质量较高。Pfam-B做为Pfam-A的补充，是一个未注释的低质量数据库，一般是由ADDA数据中的非冗余cluster自动生成的。虽然质量较低，但对于鉴定Pfam-A无法覆盖到的功能保守区域是非常有用的。
过滤 (E-value < 1e-5)
获取对应gtf文件：
直接过滤 fastq 文件即可
grep -v 非匹配项 -f 从文件读取patterns
过滤ID：
提取gtf：
提取fasta：
+ps：这里回到gtf文件提取外显子再写成fa文件
一开始个人感觉没有必要
前面已经获得了
这两个文件
在filter3_by_noncoding_exon.fa基础上根据coding.ID直接过滤就好了
（其实可以发现两种方法（fa - fa / gtf - gtf_exons - fa）得到的结果是一样的，都是142条序列）
后面发现gtf文件在进行定量以便过滤掉低表达量的lncRNA是必要的
step5：比对到NR数据库
参考：
[史上最全的BLAST本地化教程]https://zhuanlan.zhihu.com/p/556971474
[使用Diamond将宏基因组测序数据比对到Nr数据库]https://blog.csdn.net/jileiqwd/article/details/108932732
[构建NCBI本地BLAST数据库]https://wenku.baidu.com/view/7f099944c9aedd3383c4bb4cf7ec4afe04a1b1e7.html?wkts=1691388429085&bdQuery=diamond+blastx+%E5%88%B0NR%E6%95%B0%E6%8D%AE%E5%BA%93
https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/
下载速度很慢
参考：
使用ascp加速下载NR库
nohup ascp -v -k 1 -T -l 400m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nr.gz ./ &
Taxonomy ID: 9823
从NR蛋白数据库中提取Sus scrofa的蛋白质序列
首先使用TaxonKit提取特定taxons下的所有taxid
下载taxdump.tar.gz
ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
同样的方法试试拼接路径用ascp下载，成功
nohup ascp -v -k 1 -T -l 400m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/pub/taxonomy/taxdump.tar.gz ./ &
taxonkit list --ids 9823 --indent "" > pig.taxid.txt
构建一张表，第一列是taxid，后面7列跟着门纲目科属种的名称(可做可不做)
less pig.taxid.txt | taxonkit lineage | taxonkit reformat -f "{k}\t{p}\t{c}\t{o}\t{f}\t{g}\t{s}" -F | cut -f 1,3- | sed '1i\Taxid\tKingdom\tPhylum\tClass\tOrder\tFamily\tGenu\tSpecies' > pig.taxid_Ano.txt
下载taxid的accession号
https://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.FULL.gz
nohup ascp -v -k 1 -T -l 400m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/pub/taxonomy/accession2taxid/prot.accession2taxid.FULL.gz ./ &
使用csvtk在prot.accession2taxid.gz文件中提取pig.taxid所有的accession
zcat ~/refdata/prot.accession2taxid.FULL.gz | csvtk -t grep -f taxid -P pig.taxid.txt | csvtk -t cut -f accession.version >pig.taxid.acc.txt
构建NR库索引
这里所谓的“索引”并不和我们常见的fa文件fai索引一样，需要区别，这里从文件大小看来更像是另外拷贝一份加上索引信息
# 构建NR库索引
# 方法 1：使用上面下载的nr库解压后makeblastdb构建数据库
makeblastdb -in ~/database/test/nr -dbtype prot -out nr
# 方法 2：ascp 下载
ascp -v -k 1 -T -l 200m -i ~/miniconda3/envs/lncRNA/etc/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/nr.*.tar.gz ./
我们这里使用第一个方法
makeblastdb -in ~/refdata/NR/nr -dbtype prot -out nr
花了大概7个小时
使用blastdbcmd通过accession号提取NR蛋白质序列
blastdbcmd -db ~/refdata/NR/nr -entry_batch pig.taxid.acc.txt -out - | pigz -c > nr_pig.fa.gz
pigz需要自己先装一下 mamba install pigz
明明有索引文件，怎么回事：
不对，它为什么说”No alias or index file found for nucleotide database“，核酸数据库？
参考：https://zhuanlan.zhihu.com/p/556971474
查看帮助文档 blastdbcmd -help，添加 -dbtype prot 参数，还是不行：
BLAST FTP Site - BLAST® Help - NCBI Bookshelf https://www.ncbi.nlm.nih.gov/books/NBK62345/
尝试修改TITLE、移动文件目录、重新构建索引添加-parse_seqids 参数，测试：
测试用blastdbcmd从blast的数据库中抽提指定序列：KJX92028.1
参考：https://www.cnblogs.com/0820LL/p/10069841.html
blastdbcmd -db ~/refdata/NR/nr -dbtype prot -entry KJX92028.1
均不行，遂跳过：
注意这里只是一个可选过滤步骤，可以跳过，有的是不行的
我检索到的资料大部分都是我这个流程，如果大家根据这个流程走下来可以比对成功，欢迎在评论区留言
step6：过滤掉低表达量的lncRNA
通过count数量或FPKM过滤掉低表达量的lncRNA
featureCounts 统计count：
这里需要注意一个featureCounts2.0.3版本的问题：
双端测序定量添加 --countReadPairs参数
并且区别于我们之前对基因定量使用注释的参数：
nohup featureCounts -T 5 -p --countReadPairs -t exon -g gene_id -a $gtf
-o all.id.txt *bam 1>counts.id.log 2>&1 &
-t exon -g gene_id
我们这里对转录本注释定量：
-t transcript -g transcript_id
nohup featureCounts -T 12 -a ../7.lncRNA_pred/filter4_by_pfam.gtf -o raw_count.txt -p --countReadPairs -B -C -f -t transcript -g transcript_id ../3.hisat2_bams/*.bam > transcript_featureCount.log 2>&1 &
assign率是计算已经比对的reads中映射到特定特征的比例，而终端比对率是计算整个数据集中成功比对参考基因组的reads的比例
我们这里的特定特征就是那142个预测的lncRNA
# R 语言计算FPKM，筛选：FPKM > 0 in at least one sample ，得到lncRNA_id.txt
rm(list=ls())
# make count table
raw_df <- read.table(file = "raw_count.txt",header = T,skip = 1,sep = "\t")
count_df <- raw_df[ ,c(7:ncol(raw_df))]
metadata <- raw_df[ ,1:6] # 提取基因信息count数据前的几列
rownames(count_df) <- as.character(raw_df[,1])
colnames(count_df) <- paste0("SRR18059",29:64)
# calculate FPKM
countToFpkm <- function(counts, effLen)
{
N <- colSums(counts)
exp( log(counts) + log(1e9) - log(effLen) - log(N) )
}
options(scipen = 200) # 表示在200个数字以内都不使用科学计数法
fpkm = countToFpkm(count_df, metadata$Length)
# View(fpkm)
# FPKM > 0 in at least one sample
count_df.filter <- count_df[rowSums(fpkm)>0,]
write.table(rownames(count_df.filter),file="./filter6_by_fpkm_id", sep="\t",quote=F,row.names = F,col.names = F)
write.table(count_df.filter,file="./filter6_by_fpkm_mat", sep="\t",quote=F)
linux 里提取最终lncRNA的gtf和fa文件：
featureCounts对组装出的lncRNA、mRNA、other_RNA定量
我们这里主要进行lncRNA的鉴定，所以不进行其它RNA定量操作
贴上参考代码：
############## lncRNA featureaCounts ######################
nohup featureCounts -t transcript -g transcript_id \
-Q 10 --primary -s 0 -p -f -T 8 \
-a ../lncRNA.gtf \
-o ./raw_count.txt \
~/lncRNA_project/04.mapping/*.bam \
> featureCounts.log 2>&1 &
#########----------- protein_coding-------------------- #############
less -S gencode.v37.annotation.gtf | grep -w 'gene_type "protein_coding"' > protein_coding_gene.gtf
less -S gencode.v37.annotation.gtf | grep -w 'gene_type "lncRNA"' > known_lncRNA_gene.gtf
# featureCounts 定量
nohup featureCounts -t gene -g gene_id \
-Q 10 --primary -s 0 -p -f -T 8 \
-a ../protein_coding_gene.gtf \
-o ./raw_count.txt \
~/lncRNA_project/04.mapping/*.bam \
> featureCounts.log 2>&1 &
#########----------- other_RNA-------------------- #############
less -S ~/lncRNA_project/10.mRNA/gencode.v37.annotation.gtf | grep -w -v 'gene_type "protein_coding"' > other_coding_gene.gtf
# featureCounts 定量
nohup featureCounts -t gene -g gene_id \
-Q 10 --primary -s 0 -p -f -T 8 \
-a ../other_coding_gene.gtf \
-o ./raw_count.txt \
~/lncRNA_project/04.mapping/*.bam \
> featureCounts.log 2>&1 &
NONCODE v6_human数据库
blastn lncRNA.fa到 NONCODE v6_human，区分组装出的lncRNA为：know_lncRNA、novel_lncRNA
http://www.noncode.org/
NONCODE (current version v6.0) is an integrated knowledge database dedicated to non-coding RNAs (excluding tRNAs and rRNAs).
http://www.noncode.org/datadownload/NONCODEv6_human.fa.gz
nohup axel -n20 http://www.noncode.org/datadownload/NONCODEv6_human.fa.gz
官方给了md5值，下载后检查一下完整性，我现在就怀疑当时NR库下载太大了，下载的结果可能有问题，所以导致：Error: [blastdbcmd] DB contains no accession info.
删除后重新下载仍显示不完整
更换下载方式：
nohup wget -c http://www.noncode.org/datadownload/NONCODEv6_human.fa.gz &
查看下载下来的文件MD5
自己的电脑下载：
怀疑是官网MD5值没更新，继续
mamba install seqkit
Seqkit是一种用于处理和操纵序列数据的生物信息学工具。它提供了广泛的功能，如序列格式转换、修剪、过滤、统计等。
rmdup表示去除重复的序列，-s表示保留每个序列的第一个实例，-o表示将处理后的数据保存到 clean.fa文件中，-d表示将重复的序列保存到 duplicated.fa文件中，而 -D表示将关于重复序列的详细信息保存到 duplicated.detail.txt文件中
$cat ~/refdata/NONCODEv6_human.fa | seqkit rmdup -s -o clean.fa -d ddplicated.fa -D duplicated.detail.txt
[INFO] 1770 duplicated records removed
建索引
nohup makeblastdb -in clean.fa -dbtype nucl -parse_seqids -out NONCODEv6_human &
将blastn结果中e-value<=1e-10,min-identity=80%,min-coverage=50%的序列筛选出来当作比对上的序列
nohup blastn -db NONCODEv6_human -evalue 1e-10 -num_threads 10 -max_target_seqs 5 -query ../8.featureCounts/lncRNA.fa -outfmt ' 6 qseqid sseqid pident qcovs length mismatch gapopen qstart qend sstart send qseq sseq evalue bitscore' -out lncRNA.txt
（这里比对NONCODE成功了，以后如果还是要做比对NR，可以在保证文件完整性的前提下，根据这里比对NONCODE一些步骤来debug）
过滤这些匹配上的id，提取余下lncRNA的gtf和fa：
根据匹配上NONCODE已知lncRNA的过滤掉6个
lncRNA的功能推断
大量lncRNA的功能是未知的，但是它们主要是cis-regulators，所以可以根据它们
临近的蛋白编码基因功能来近似推断 然后表达量的相关性也可以类推到
根据位置关系推断
nohup bedtools window -a ../9.NONCODE_blastn/lncRNA.gtf -b ~/refdata/gencode.v44.annotation_RmChr.gtf -l 100000 -r 10000 > lncRNA_locus.txt &
输入文件的每个区域周围创建一个窗口，左侧扩展100,000个单位，右侧扩展10,000个单位，并将结果保存
locus：
lncRNA：
表达量的相关性
以上，就是本期全部内容
https://mp.weixin.qq.com/s/xqO7Wur2cMLLwRUSgR1yuQ