lncRNA数据处理（get fasta>>>novel lncRNA))）

情况说明：我有两个实验，每个实验两个重复，我现在要鉴定到新的lncRNA之后进行差异分析。大体思路，我会先通过常规的鉴定novel lncRNA方法鉴定到新的转录本，将转录本与参考基因组的转录本合并，之后进行定量分析。因此，我之前处理好各个文件的接头序列之后比对基因组后，会将四个bam文件合并成一个，之后进行鉴定（通过length、fpkm、cov、编码能力预测）之后会将鉴定到的转录本合并（与基因组的转录本），使用RSEM定量

使用到的软件有--fastQC、fastp、STAR、stringtie、cuffcompare、cpc（cpc2）、CNCI、CPAT、

1.常规的数据处理步骤

1）拿到数据后先进行质检（常用fastqc看测序质量、用fastp做处理去接头什么的）

fastqc -t 4 ` ls * `

fastp -i ${file}-r1.fastq.gz -o ${file}-f.r1.fastq.gz -I ${file}-r2.fastq.gz -O ${file}-f.r2.fastq.gz -Q --thread=5 -5 --detect_adapter_for_pe 2> ${file}.log

###${file}是我用for循环做的变量（这里应该是你的测序文件，PE150双端150bp测序，所以会有两个文件，分别是r1和r2。fastp软件使用可以直接看manual ，比较杂，我主要用了其中的去接头处理，关闭质量过滤”Q” 2>log文件是使用该软件处理过程中生成的日志。

2）开始进行比对（根据物种来选择一个合适的比对软件，植物多倍体推荐bwa，STAR，我用的是STAR，大体参数有这些

STAR --runThreadN 20 \

# --outSAMmapqUnique 60 \

# --twopassMode Basic \

# --outSAMtype BAM SortedByCoordinate \

# --limitBAMsortRAM 200000000000 \

# --outFilterMultimapNmax 1 \

# --outSAMprimaryFlag AllBestScore \

# --outReadsUnmapped Fastx \

# --outFilterMatchNmin 50 \

# --alignMatesGapMax 10000 \

# --outSAMstrandField intronMotif \

# --outFilterMismatchNoverLmax 0.03 \

# --outFilterType BySJout \

# --readFilesCommand zcat \

# --alignSoftClipAtReferenceEnds No \

# --alignEndsType Local \

# --alignIntronMin 20 \

# --alignIntronMax 15000 \

# --outSAMattributes All \

# --outFilterMismatchNoverReadLmax 0.02 \

这个参数比较多，但实际应用的时候可以相应减少，其实很多时候默认参数就可以，发文章的时候默认参数比较有说服力。

--runThreadN 20 --genomeDir $index \

--readFilesIn ${name}_1P.fq.gz ${name}_2P.fq.gz \

--outReadsUnmapped Fastx \

--readFilesCommand zcat \

--outSAMtype BAM SortedByCoordinate \

--outFileNamePrefix ./${name}_star/${name}

###这是我用的参数，仅供参考，我是为了之后找lncRNA

3)比对完之后要进行组装，有cufflinks stringtie等，此处用的是stringtie

"""stringtie -p 4 -G ref-genome.gtf -o output.gtf input_file """

###有的是在这里会将多个组装好的gtf文件merge到一起，可以使用stringtie 的merge

4）使用cuffconpare进行分类

""" cuffcompare -r ref-genome.gtf -o output input_file """

###会从中得到好多文件，一般我们会从*.tmap文件中进行筛选

"""awk '{if($7>=0.5 && $10 > 1 && $11 >200) print $0}' 自己命名的.tmap |awk '{if ($3=="u" || $3=="x" || $3=="i" || $3=="j" || $3=="o"){print $0}}' > newfile.gtf

###其中$7是fpkm值、$10是cov值、$11是length ，使用管道直接筛选出$3是u、x、i、j、o的code，具体的分类值可以从cuffcompare的manual看到

5）到了写脚本的时候啦，要做的有，从上一步筛选到的转录本的到它的转录本序列（我是得到它的每个转录本有多少个exon，位置是多少，从已知基因组中得到序列后将其合并到一起，形成转录本序列，之后进行coding potential calculator（CPC）CNCI等软件来做，CPC我是安装到liunx 本地的，