前面我们介绍了一系列的LncRNA鉴定相关文献的案例精选：

- 4个发育时间点的总共12个鸡转录组测序样本的长非编码RNA的鉴定

• 59匹马的8个组织的长非编码RNA的鉴定
• 9个组织的37个样本的大豆的长非编码RNA的鉴定
• 中国荷斯坦奶牛新的lncRNA全基因组鉴定
• 纯数据挖掘之仔猪的长非编码RNA的鉴定
• 高原牦牛的长非编码RNA的鉴定
• 糖尿病视网膜病变患者的长非编码RNA的鉴定
• 苹果小卷蛾长非编码RNA的鉴定
• 胡萝卜的长非编码RNA的鉴定
• 急性髓系白血病的lncRNAs表观遗传图谱

相信这个时候大家对LncRNA鉴定的整体流程应该是有所了解了，这里就不再赘述，很多小伙伴留言找这个代码，现在把上下游分开整理给大家。

这里说的上游分析主要是从fastq测序文件开始，走hisat2和stringtie流程

raw 数据简介

保存在 /home/data/ccrcc/public_data/ccrcc/目录下，节选 如下所示:

4.0G 1月 6 09:49 ccrcc/SRR10744251_1.fastq.gz
4.4G 1月 6 09:52 ccrcc/SRR10744251_2.fastq.gz
4.2G 1月 6 09:55 ccrcc/SRR10744252_1.fastq.gz
4.6G 1月 6 10:03 ccrcc/SRR10744252_2.fastq.gz

fastqc质控

我使用的批量命令是：

ls $HOME/lncRNA_project/01.raw_data/*.gz > config

cat > 02.fastqc.sh
config=$1
number1=$2
number2=$3
cat $1 | while read id
do
 if((i%$number1==$number2))
 then
 fastqc -t 5 $id -o ./ 
 fi 
 i=$((i+1))
done

for i in {0..4}
do
(nohup bash 02.fastqc.sh config 5 $i 1>log.$i.txt 2>&1 & )
done

得到 $HOME/lncRNA_project/02.fastqc/01.raw_qc/ ,目录如下文件：

592K 2月 5 18:53 SRR10744251_1_val_1_fastqc.html
322K 2月 5 18:53 SRR10744251_1_val_1_fastqc.zip
594K 2月 5 18:54 SRR10744251_2_val_2_fastqc.html
322K 2月 5 18:54 SRR10744251_2_val_2_fastqc.zip

接下来对raw的fastq文件进行trim

我使用的命令是：

ls $HOME/lncRNA_project/01.raw_data/*_1.fastq.gz > 1
ls $HOME/lncRNA_project/01.raw_data/*_2.fastq.gz > 2
paste 1 2 > config

cat > 03.trim.sh
config=$1
number1=$2
number2=$3
cat $1 | while read id
do
 if((i%$number1==$number2))
 then
 arr=(${id})
 fq1=${arr[0]}
 fq2=${arr[1]}
 trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o ./ $fq1 $fq2
 fi ## end for number1
 i=$((i+1))
done

for i in {0..9}
do
(nohup bash 03.trim.sh config 10 $i 1>log.$i.txt 2>&1 & )
done

$HOME/lncRNA_project/ 得到的文件如下:

 3.9G 2月 4 14:09 03.trim/SRR10744251_1_val_1.fq.gz
 4.3G 2月 4 14:09 03.trim/SRR10744251_2_val_2.fq.gz
 4.1G 2月 4 14:12 03.trim/SRR10744252_1_val_1.fq.gz
 4.4G 2月 4 14:12 03.trim/SRR10744252_2_val_2.fq.gz

hisat2比对

我使用的命令是：

mkdir 04.mapping
ls $HOME/lncRNA_project/03.trim/*_1.fq.gz > 1
ls $HOME/lncRNA_project/03.trim/*_2.fq.gz > 2
ls $HOME/lncRNA_project/03.trim/*_1.fq.gz | cut -d "/" -f 7 | cut -d "_" -f 1 > 0
paste 0 1 2 > config

cd 04.mapping

cat > 04.hg38_mapping.sh
index=/home/data/server/reference/index/hisat/hg38/genome
config=$1
number1=$2
number2=$3
cat $1 | while read id
do
 if((i%$number1==$number2))
 then
 arr=(${id})
 fq1=${arr[1]}
 fq2=${arr[2]}
 sample=${arr[0]}
 hisat2 -p 2 -x ${index} -1 $fq1 -2 $fq2 | samtools sort -@ 6 -o ./${sample}.bam -
 fi 
 i=$((i+1))
done

for i in {0..6}
do
(nohup bash 04.hg38_mapping.sh config 7 $i 1>log.$i.txt 2>&1 & )
done

得到的文件节选如下所示：

6.1G 2月 6 14:13 04.mapping/SRR10744251.bam
6.5G 2月 6 14:12 04.mapping/SRR10744252.bam
8.1G 2月 6 14:50 04.mapping/SRR10744253.bam
6.6G 2月 6 14:38 04.mapping/SRR10744254.bam

samtools index

前面仅仅是比对，把fastq文件转为了bam文件，并没有对bam文件构建index，后续很多操作bam文件的命令和软件都依赖于bam文件的index文件。

我使用的代码：

ls $HOME/lncRNA_project/04.mapping/*.bam > 1
ls $HOME/lncRNA_project/04.mapping/*.bam | cut -d "/" -f 7 | cut -d "_" -f 1 | cut -d "." -f 1 > 0
paste 0 1 > config

cat > index.sh

config=$1
number1=$2
number2=$3
cat $1 | while read id
do
 if((i%$number1==$number2))
 then
 arr=(${id})
 input=${arr[1]}
 sample=${arr[0]}
 samtools index -@ 4 ${input} ./${sample}.bam.bai
 fi ## end for number1
 i=$((i+1))
done

for i in {0..4}
do
(nohup bash index.sh config 5 $i 1>log.$i.txt 2>&1 & )
done

得到如下文件，节选：

4.8M 2月 8 20:33 04.mapping/SRR10744251.bam.bai
5.1M 2月 8 20:33 04.mapping/SRR10744252.bam.bai
5.9M 2月 8 20:33 04.mapping/SRR10744253.bam.bai
5.1M 2月 8 20:33 04.mapping/SRR10744254.bam.bai

bam 转 bw

因为前面的bam文件过于巨大，不方便载入igv查看，所以bam 转 bw，然后可以下载到本地电脑载入igv啦。

使用脚本

ls $HOME/lncRNA_project/04.mapping/*.bam > 1
ls $HOME/lncRNA_project/04.mapping/*.bam | cut -d "/" -f 7 | cut -d "_" -f 1 | cut -d "." -f 1 > 0
paste 0 1 > config

cat > bw.sh

config=$1
number1=$2
number2=$3
cat $1 | while read id
do
 if((i%$number1==$number2))
 then
 arr=(${id})
 input=${arr[1]}
 sample=${arr[0]}
 bamCoverage -p 2 -b ${input} -o ./${sample}.bw
 fi ## end for number1
 i=$((i+1))
done

for i in {0..4}
do
(nohup bash bw.sh config 5 $i 1>log.$i.txt 2>&1 & )
done

得到如下文件：

 87M 2月 8 20:49 04.mapping/bw/SRR10744251.bw
 95M 2月 8 20:51 04.mapping/bw/SRR10744252.bw
 99M 2月 8 20:53 04.mapping/bw/SRR10744253.bw
 79M 2月 8 20:52 04.mapping/bw/SRR10744254.bw

使用stringtie软件的assmebly功能

使用的命令：

ls $HOME/lncRNA_project/04.mapping/*.bam > 1
ls $HOME/lncRNA_project/04.mapping/*.bam | cut -d "/" -f 7 | cut -d "_" -f 1 | cut -d "." -f 1 > 0
paste 0 1 > config

cat > 05.stringtie.sh

config=$1
number1=$2
number2=$3
gtf=$HOME/reference/human/gtf/gencode.v37.annotation.gtf
cat $1 | while read id
do
 if((i%$number1==$number2))
 then
 arr=(${id})
 input=${arr[1]}
 sample=${arr[0]}
 stringtie -B -p 4 -G $gtf -o ./${sample}.stringtie.gtf -A ./${sample}.tab -l ${sample} $input
 fi ## end for number1
 i=$((i+1))
done
# 需要理解 -B 参数哦 ， 如果StringTie与-B选项一起运行,它将返回Ballgown输入文件

for i in {0..4}
do
(nohup bash 05.stringtie.sh config 5 $i 1>log.$i.txt 2>&1 & )
done

得到如下文件：

200M 2月 17 14:06 05.stringtie/01.stringtie_gtf/SRR10744251.stringtie.gtf
198M 2月 17 14:07 05.stringtie/01.stringtie_gtf/SRR10744252.stringtie.gtf
208M 2月 17 14:09 05.stringtie/01.stringtie_gtf/SRR10744253.stringtie.gtf
197M 2月 17 14:07 05.stringtie/01.stringtie_gtf/SRR10744254.stringtie.gtf

merge gtf文件

因为每个样品都输出了自己的gtf文件，后面进行de novo 的lncRNA鉴定的时候其实是合并全部的样品，所以这里的gtf文件也需要合并。

我是用的代码是：

ls ../01.stringtie_gtf/*.gtf > mergelist.txt
gtf=$HOME/reference/human/gtf/gencode.v37.annotation.gtf
nohup stringtie --merge -p 8 -G $gtf -o ./stringtie_merged.gtf ./mergelist.txt > merge.log 2>&1 &

得到如下的文件：

8.5K 2月 17 20:23 mergelist.txt
 22 2月 17 20:26 merge.log
536M 2月 17 20:44 stringtie_merged.gtf

Stringtie估算转录本表达量

这个时候可以使用我们构建好了全新的gtf文件，对全部的bam文件进行定量。

仅仅是跑个程序而已，其实这个表达量文件，我们后续也可以不使用，仍然是传统的featureCounts定量即可。

ls $HOME/lncRNA_project/04.mapping/*.bam > 1
ls $HOME/lncRNA_project/04.mapping/*.bam | cut -d "/" -f 7 | cut -d "_" -f 1 | cut -d "." -f 1 > 0
paste 0 1 > config

cat > 05.estimate_abundance

config=$1
number1=$2
number2=$3
cat $1 | while read id
do
 if((i%$number1==$number2))
 then
 arr=(${id})
 input=${arr[1]}
 sample=${arr[0]}
 stringtie -e -B -p 4 -G ~/lncRNA_project/05.stringtie/02.merge_gtf/stringtie_merged.gtf -o ./ballgown/${sample}/${sample}.gtf $input

fi ## end for number1
 i=$((i+1))
done

for i in {0..4}
do
(nohup bash 05.estimate_abundance config 5 $i 1>log.$i.txt 2>&1 & )
done

这个时候的hisat2和stringtie流程就结束了，但是组装好的gtf文件还需要进行一系列的评估和过滤。