鉴定新的lncRNA之上游流程

Posted on by

好奇怪哦，我们前面的 lncRNA-seq数据分析之新lncRNA鉴定和注释视频课程众筹 ，感兴趣的人似乎不多额，免费的啊，既然感兴趣人不多，这个视频课程就取消免费了哈！那个群大家仍然是可以进入，拿资料和代码，视频我就不录制了，感兴趣的人不多，我搞自媒体图的就是一个名声，没有人看，我浪费时间干嘛！

上游流程，通常指的是ngs测序数据fastq文件，在服务器级别的计算资源里面的一系列处理。因为个人电脑很难hold住，而且流程很少变动，所以通常是公司代替客户完成。属于吃力不讨好的技能，学习成本极高，但是学完之后用的频率很低！下游分析流程，指的是拿到上游分析结果之后在自己个人电脑里面的统计可视化图表制作。

比如RNA-seq数据，上游就是fastq的质量控制，比对，定量，最后拿到表达矩阵。而下游就是表达矩阵的一系列统计学分析， 包括PCA，相关性热图，层次聚类图，差异分析，火山图，表达量热图，GO/KEGG数据库功能注释等等。

对lncRNA-seq鉴定新的lncRNA流程来说，也是上下游独立。我们先介绍一下上游流程：

首先使用conda来创建LncRNA-seq的实战软件环境

conda create -n lncRNA 
conda activate lncRNA 
conda install -y -c bioconda hisat2 stringtie samtools fastp gffcompare
# conda search gffcompare

然后下载猪的参考基因组fasta序列，并且构建hisat2的索引文件

参考：猪狗的参考基因组构建索引

mkdir -p ~/reference/genome/
cd ~/reference/genome/
mkdir pig
cd pig 
wget ftp://ftp.ensembl.org/pub/release-99/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz 
gunzip -d Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz 
# 需要注意文件大小，以及参考基因组是否下载成功哦！

wget ftp://ftp.ensembl.org/pub/release-99/gtf/sus_scrofa/Sus_scrofa.Sscrofa11.1.99.chr.gtf.gz
gunzip -d Sus_scrofa.Sscrofa11.1.99.chr.gtf.gz

nohup hisat2-build -p 4 Sus_scrofa.Sscrofa11.1.dna.toplevel.fa pig_hisat2 &

一定要注意 hisat2-build 根据参考基因组的fasta序列文件构建索引对内存的消耗，以及输出文件的大小哦。

接着下载文章lncRNA-seq的fastq测序数据

参考：使用ebi数据库直接下载fastq测序数据

cat fq.txt |while read id
do
ascp -QT -l 300m -P33001 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh \
era-fasp@$id .
done
# nohup bash step1-aspera.sh 1>step1-aspera.log 2>&1 &

这个脚本会根据你在EBI里面搜索到的 fq.txt 路径文件，来批量下载fastq测序数据文件。这个 fq.txt 文件里面的路径还是蛮有规律的，如下：

fasp.sra.ebi.ac.uk:/vol1/fastq/SRR180/009/SRR1805929/SRR1805929_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR180/009/SRR1805929/SRR1805929_2.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR180/000/SRR1805930/SRR1805930_1.fastq.gz
fasp.sra.ebi.ac.uk:/vol1/fastq/SRR180/000/SRR1805930/SRR1805930_2.fastq.gz

值得注意的是，EBI数据库也不是那么稳定，有时候这个下载会失败，大家可以等候几天重新尝试，或者干脆转到去sra数据库下载。

批量运行fastp+hisat2+stringtie流程

需要学习和使用3个软件的用法，fastp+hisat2+stringtie，理解参数：

conda activate lncRNA
index=/home/yb77613/reference/genome/pig/pig_hisat2
gtf=/home/yb77613/reference/genome/pig/Sus_scrofa.Sscrofa11.1.99.chr.gtf

fastp -i 1.raw_fq/${id}_1.fastq.gz \
 -o 2.clean_fq/${id}_1.fastp.fq.gz \
 -I 1.raw_fq/${id}_2.fastq.gz \
 -O 2.clean_fq/${id}_2.fastp.fq.gz \
 -l 36 -q 20 --compression=6 \
 -R ${id} -h ${id}.html
fq1=2.clean_fq/${id}_1.fastp.fq.gz
fq2=2.clean_fq/${id}_2.fastp.fq.gz
hisat2 -p 4 -x $index -1 $fq1 -2 $fq2 | \
samtools sort -@ 4 -o 3.hisat2_bams/$sample.bam -

stringtie -p 4 -G $gtf \
 -o 4.stringtie_gtfs/$sample.gtf \
 -l $sample 3.hisat2_bams/$sample.bam

批量处理！

合并gtf后鉴定和筛选新的lncRNA

conda activate lncRNA 
ls `pwd`/4.stringtie_gtfs/*gtf > 5.lncRNA/gtf.list
cd 5.lncRNA

gtf=/home/yb77613/reference/genome/pig/Sus_scrofa.Sscrofa11.1.99.chr.gtf
nohup stringtie --merge -p 6 -G $gtf -o stringtie_merged.gtf gtf.list > merge.log 2>&1 &

gff=Sus_scrofa.Sscrofa11.1.93.gff3
gffcompare -r $gff -o gffcomp93 stringtie_merged.gtf
# 49448 reference transcripts loaded.
# 69 duplicate reference transcripts discarded.
# 113472 query transfrags loaded.

gtf=/home/yb77613/reference/genome/pig/Sus_scrofa.Sscrofa11.1.99.chr.gtf
gff=$gtf
gffcompare -r $gff -o gffcomp99 stringtie_merged.gtf
# 60954 reference transcripts loaded.
# 84 duplicate reference transcripts discarded.
# 113472 query transfrags loaded.

可以看到，使用不同版本的gtf文件，对结果的影响很明显哦，鉴定新的lncRNA，也许2016年你认为是新的lncRNA，到2018年就已经被记录在册啦。

如果是人类数据，我们通常是使用gencode的gtf文件哦

gtf=/home/yb77613/reference/gtf/gencode/gencode.v32.annotation.gtf
nohup stringtie --merge -p 6 -G $gtf -o stringtie_merged.gtf gtf.list > merge.log 2>&1 &
gffcompare -r $gtf -o gffcomp stringtie_merged.gtf

其它物种，大家类推即可，主要是人类和小鼠，做这个流程的意义不是很大。当然，如果你强行跑流程，也是可以的，结果的生物学解释就比较考验人。

我们会有一系列lncRNA-seq数据分析相关资料共享，前面的 lncRNA-seq数据分析之新lncRNA鉴定和注释视频课程众筹 有入群交流方式哈。（因为感兴趣人不多，所以我就不录制视频啦，浪费时间，还赚不到名声。）

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

​

Comments are closed.