LncRNA鉴定

0.简介

主要参考一篇NC文章鉴定lncRNA流程。

https://www.nature.com/articles/s41467-018-07500-7

1.质控

#!/bin/bash
for i in `ls *.fastq`;do x=${i/.fastq/};


echo $x
##Trimmomatic
java -jar /biosoft/Trimmomatic-0.38/trimmomatic-0.38.jar SE -threads 4 -phred33 $x.fastq $x.1.fastq ILLUMINACLIP:/biosoft/Trimmomatic-0.38/adapters/TruSeq2-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
##更改adapter文件为TruSeq2-SE.fa，使得fastqc Adapter Contenr 项达标

##
musket -k 21 50000000 $x.1.fastq -o $x.2.fastq

fastqc $x.2.fastq
done

2.hisat2+stringtie

使用hista2+stringtie速度明显快于tophat+cufflinks。

hisat2-build 构建参考序列索引

 hisat2-build ~/lncRNA/Arabidopsis_thaliana.TAIR10.fa genome

hisat2+stringtie流程：

#!/bin/bash
for i in `ls *.2.fastq`;do s=${i/.2.fastq/};
echo $s

##hisat2
hisat2 -x /home/manager/lncRNA/genome -U $s.2.fastq --dta -S $s.sam

##samtools 将.sam转为.bam
samtools view -bS $s.sam >$s.bam

##samtools sort 排序
samtools sort $s.bam $s.sort

##tringtie
stringtie $s.sort.bam --rf -l $s -G ../Arabidopsis_thaliana.TAIR10.47.gff3 -o $s.gtf

done

tophat+cufflinks（备选）：

##1.bowtie
bowtie2-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa genome

##2.tophat
tophat --library-type fr-firststrand -g 5 -G Arabidopsis_thaliana.TAIR10.47.gff3 genome 2_SRR627952.fastaq

##3.cufflinks
cufflinks -g Arabidopsis_thaliana.TAIR10.47.gff3 --min-frags-per-transfrag 1   accepted_hits.bam

整合gtf

cuffmerge -g ~/lncRNA/Arabidopsis_thaliana.TAIR10.47.gff3 -s ../Arabidopsis_thaliana.TAIR10.fa  gtf.txt

#或者stringtie--merge：
stringtie --merge -o merged.gtf -G ../Arabidopsis_thaliana.TAIR10.47.gff3  *.gtf

获得转录本丰度值:

#!/bin/bash
for i in `ls *.sort.bam`;do x=${i/.sort.bam/};
echo $x
stringtie -G merged.gtf -e -A $x.A.txt $x.sort.bam
done

整合丰度值文件:

perl integrate_fpkm.pl >fpkm.txt

生成文件：fpkm.txt

fpkm.txt是个两列文件，第一列是gene_id，第二列是fpkm值

3.筛选lncRNA

3.1.保留 class codes 为 "x"，"i"，"u"的转录本

awk '/class_code "u"/||/class_code "i"/||/class_code "x"/' merged.gtf >filter1.gtf

3.2.过滤长度小于150nt,fpkm值小于1的数据

perl filter_gtf.pl fpkm.txt filter1.gtf >filter2.gtf

根据gtf提取序列

gffread  filter.fa -g ../Arabidopsis_thaliana.TAIR10.fa filter2.gtf

3.3.去除有蛋白质编码潜能的数据

3.3.1.blastx

与swiss-prot的植物蛋白质数据库blastx比对，截取掉 e-value < 10-4，和有高匹配度的数据（alignment length ≥40 aa, percent identity ≥35%）

首先处理下载的数据库文件

cat uniprot_sprot_plants.dat |awk '{if (/^ /) {gsub(/ /, ""); print} else if (/^AC/) print ">" $2}' |sed 's/;$//'> protein.fa

blastx

##建库
makeblastdb -in protein.fa -dbtype prot -out sprot_plant
##blast
blastx -db sprot_plant -out blastx_out.txt -outfmt 6 -query ../filter.fa -num_threads 4

生成文件：blastx_out.txt

3.3.2.CNCI

原文中使用CPC，但由于CPC的官网我无法访问，因此暂由CNCI替代

python /home/manager/biosoft/CNCI-master/CNCI.py -f filter.fa -o CNCI -m pl -p 4

生成文件：CNCI.index

3.3.3.整合

整合blastx和cnci筛选得到的数据

perl Integrate_filter.pl blastx/blastx_out.txt CNCI/CNCI.index >lncRNA_id.txt

得到文件：lncRNA_id.txt 该文件中包含提出的lncRNA ID

#脚本

#1.integrate_fpkm.pl

use strict;
use warnings;
my %hash;
my @pm=glob '*.A.txt';
foreach (@pm){
open TXT,"$_" or die "erro : cannot open the file";
while(my $line = <TXT>){
	my ($gene_id,$fpkm)=(split /\s/,$line)[0,7];
	if(!exists $hash{$gene_id}){
		$hash{$gene_id}=$fpkm;	
	}else{
		if($fpkm>$hash{$gene_id}){
			$hash{$gene_id}=$fpkm;
	}
	}
}
close TXT;
}


foreach my $key (keys %hash){
	print $key."\t".$hash{$key}."\n";
}

#2.Integrate_filter.pl

use strict;
use warnings;
my %hash;
my $txt=$ARGV[0];
my $CNCI=$ARGV[1];
open TXT,"$txt"or die "cannot open the blast file";


while(my $line = <TXT>){
	chomp $line;
	my($id,$identity,$length,$e_value)=(split /\s/,$line)[0,2,3,10];
	if($e_value<1e-3){next;}
	if($length>=40 && $identity>=35){next;}
	if(!exists $hash{$id}){
		$hash{$id}=1;
	}
	
}
close TXT;


open CNCI,"$CNCI"or die "cannot open the CNCI file";
while(my $line2=<CNCI>){
	chomp $line2;
	my($id,$index)=(split /\s/,$line2)[0,1];
	if($index eq 'noncoding'){
		if(exists $hash{$id}){
			print"$id\n";
	}
	}
}
close CNCI;