LncRNA鉴定上游分析

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

相信这个时候大家对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文件还需要进行一系列的评估和过滤。

 很多小伙伴留言表示我们的课程里面的下游分析太少了，明天我们就发布全套下游分析代码哦！

赞 (0)

相关推荐

  • 4 比对到参考基因组输出bam文件

    进到align目录 对质量好的测序数据进行比对 1. 一个个比对,生成BAM文件 align目录 sample=SRR7696207 bwa mem -t 2 -R "@RG\tID:$sa ...

  • 完美 | GXF Fix 修复 / 优化基因结构注释信息文件 - GTF/GFF3

    写在前面 目前基因组测序和组装成本几乎已经到任何一个课题组都可以单独负担的价码,大量物种的基因组序列被测定和释放.与此同时,对应的基因结构注释信息文件,如GTF或GFF3文件等,也可公开下载. 对于绝 ...

  • 鉴定新的lncRNA之上游流程

    好奇怪哦,我们前面的 lncRNA-seq数据分析之新lncRNA鉴定和注释视频课程众筹 ,感兴趣的人似乎不多额,免费的啊! 既然感兴趣人不多,这个视频课程就取消免费了哈! 那个群大家仍然是可以进入, ...

  • 一作解读| 两对小麦分蘖近等基因系的lncRNA鉴定分析

    长链非编码RNA(lncRNA)作为新兴的转录后调控因子,已在拟南芥.水稻.玉米等多种作物中被广泛鉴定,可通过顺式或反式作用参与植物逆境响应.有性生殖.器官形成等多种生长发育活动.分蘖是小麦生长发育的 ...

  • 畜禽养殖废弃物资源化利用装备国推鉴定情况分析

    近年来,国家相继出台了多项政策支持发展畜禽养殖废弃物资源化利用,市场需求旺盛,产业发展迅速,产品种类和生产企业数量都不断增长,规模养殖场畜禽养殖废弃物资源化利用装备配套比例快速提高. 广大农业从业人员 ...

  • 把vcf文件转换为maf格式，肿瘤外显子上游分析教程到此为止

    完全看懂这个教程需要的背景知识比较多! 生信技能树GATK4系列教程 GATK4的gvcf流程 你以为的可能不是你以为的 新鲜出炉的GATK4培训教材全套PPT,赶快下载学习吧 曾老湿最新私已:GAT ...

  • lncRNA-seq数据分析之新lncRNA鉴定和注释视频课程众筹

    前面我系统性的总结了:lncRNA的一些基础知识 ,和lncRNA芯片的一般分析流程 ,还有LncRNA-seq的一般分析流程 ,里面提到了一个目前非常小众的分析方向,就是新lncRNA鉴定和注释,因 ...

  • 看优秀本科生如何一周内学会Linux进而搞定RNA-seq上游分析

    距离公布要带500个优秀本科生入门生物信息学的活动不到一个月,虽然真正入选不到一百,但是培养成绩喜人,出勤率接近百分之百,大部分人在短短两个星期就完成了R基础知识学习,Linux认知,甚至看完了转录组 ...

  • 明码标价之普通转录组上游分析

    最近有粉丝在我们<生信技能树>公众号后台付费求助,想重新分析一下他自己领域的一个文章的转录组测序数据.因为那个文章并没有提供表达量矩阵,不过万幸的是人家上传了原始测序数据. 因为他的课题是 ...

  • 单细胞转录组上游分析之shell回顾

    课程笔记 目录 单细胞转录组学习笔记-1 单细胞转录组学习笔记-2 前言 目前主流的单细胞测序技术主要有两种:主打基因数量的smart-seq2和主打细胞数量的10X Genomics.单细胞转录组分 ...

  • 单细胞转录组数据处理之上游分析流程

    如果你想亲手分析自己的数据,生信技能树联合单细胞天地也推出了:7个小时的单细胞转录组视频课程(限时免费)  足足7个小时,34集,涵盖了单细胞转录组背景知识以及数据处理的知识体系,希望能够帮助到你的课 ...