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

Posted 2023-04-08

参考技术A 情况说明：我有两个实验，每个实验两个重复，我现在要鉴定到新的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 本地的，

python处理fasta文件，ID和序列放在一行

#!/usr/bin/python
#-*- coding:utf-8 -*-
"处理fasta文件，将ID号和序列放在一行"
import sys
with open(sys.argv[1]) as f:
    fw=open(‘out.fasta‘, ‘w‘)
    line=f.read()
    line=line.replace(‘\n‘, ‘‘).replace(‘>‘, ‘\n>‘)
    for aa in line:
        fw.write(aa)
    fw.close()
"""
>chr1|hos107.1#gene1
ACACTCCCGGGCCCCCCCCCCCC
ACCTTTCAAAAAAAAAAAAAAA
AATTTTCCCCCCAAAGGGG
>chr1|hos107.2#gene2
ACACTCCCGGGCCCCCCCCCCCC
ACCTTTCAAAAAAAAAAAAAAA
AATTTTC
>chr1|hos107.4#gene3
ACACTCCCGGGCCCCCCCCCCCC
ACCTTTCAAAAAAAAAAAAAAA
AATTTTC
>chr1|hos107.5#gene4
ACACTCCCGGGCCCCCCCCCCCC
ACCTTTCAAAAAAAAAAAAAAA
AATTTTC
"""
"""
>chr1|hos107.1#gene1ACACTCCCGGGCCCCCCCCCCCCACCTTTCAAAAAAAAAAAAAAAAATTTTCCCCCCAAAGGGG
>chr1|hos107.2#gene2ACACTCCCGGGCCCCCCCCCCCCACCTTTCAAAAAAAAAAAAAAAAATTTTC
>chr1|hos107.4#gene3ACACTCCCGGGCCCCCCCCCCCCACCTTTCAAAAAAAAAAAAAAAAATTTTC
>chr1|hos107.5#gene4ACACTCCCGGGCCCCCCCCCCCCACCTTTCAAAAAAAAAAAAAAAAATTTTC
"""

#提取目标序列
f=open(‘./out.fasta‘, ‘r‘)
fw=open(‘target.fasta‘, ‘w‘) 
for line in f.readlines():
    if line.startswith(‘>chr1|hos107.1‘):
        fw.write(line)
f.close()
fw.close()


"""可以从上述处理好的单行文件out.fasta中提取指定目标ID的文件，并将其
写入到target.fasta文件中"""

#整体思路：
#先统一fasta文件格式从test.fasta----out.fasta
#取出目标ID序列：out.fasta----target.fasta