2023-07-21使用Trinity进行转录组De novo组

作者: 麦冬花儿 | 来源:发表于2023-08-12 09:57 被阅读0次

转录组入门学习（四）
使用Trinity 软件进行拼接（无参转录组）
全基因组测序从头测序(de novo sequencing)
转录组数据库的基本使用（一）-GO数据库
转录组组装软件：Trinity (转载)
空间转录组分析-工具箱
转录组测序技术和结果解读（五）——参考基因组mapping
基因组从头组装软件
基因组结构注释软件列表
4、RNAseq(4)--Hisat2进行序列比对及Samtoo

软件安装与下载

# Installing trinity (https://github.com/trinityrnaseq/trinityrnaseq/releases)
#wget https://github.com/trinityrnaseq/trinityrnaseq/releases/download/v2.15.1/trinityrnaseq-v2.15.1.FULL.tar.gz -P ~/software/
tar zxf ~/software/trinityrnaseq-v2.15.1.FULL.tar.gz
mv trinityrnaseq-v2.15.1/ /opt/biosoft/Trinity-v2.15.1
cd /opt/biosoft/Trinity-v2.15.1
make -j 4
make plugins
echo 'PATH=$PATH:/opt/biosoft/Trinity-v2.15.1/' >> ~/.bashrc
source ~/.bashrc

# Installing Jellyfish (http://www.genome.umd.edu/jellyfish.html)
# Trinity的运行需要依赖jellyfish version 2。Jellyfish两个版本并存，注意不要下错了。
#wget https://github.com/gmarcais/Jellyfish/releases/download/v2.3.0/jellyfish-linux -P ~/software/
cp /home/train/software/jellyfish-linux /opt/biosoft/Trinity-v2.15.1/jellyfish
chmod 755 /opt/biosoft/Trinity-v2.15.1/jellyfish

# Installing samlomon (https://combine-lab.github.io/salmon/ | https://github.com/COMBINE-lab/salmon)
#wget https://github.com/COMBINE-lab/salmon/releases/download/v1.10.0/salmon-1.10.0_linux_x86_64.tar.gz -P ~/software/
tar zxf ~/software/salmon-1.10.0_linux_x86_64.tar.gz -C /opt/biosoft/
echo 'PATH=$PATH:/opt/biosoft/salmon-latest_linux_x86_64/bin/' >> ~/.bashrc
source ~/.bashrc

#tar zxf /home/train/software/R-4.3.1.Rocky9.2_opt_sysoft.tar.gz -C /opt/sysoft/
#echo 'PATH=$PATH:/opt/sysoft/R-4.3.1/bin/' >> ~/.bashrc
#source ~/.bashrc

#Installing RSEM (http://deweylab.github.io/RSEM/)
#wget https://github.com/deweylab/RSEM/archive/v1.3.3.tar.gz -O ~/software/RSEM-v1.3.3.tar.gz
tar zxf ~/software/RSEM-v1.3.3.tar.gz -C /opt/biosoft/
cd /opt/biosoft/RSEM-1.3.3/
make -j 4
echo 'PATH=$PATH:/opt/biosoft/RSEM-1.3.3/' >> ~/.bashrc
source ~/.bashrc 

#Installing kallisto (https://github.com/pachterlab/kallisto/releases)
#wget https://github.com/pachterlab/kallisto/releases/download/v0.50.0/kallisto_linux-v0.50.0.tar.gz -P ~/software/
tar zxf ~/software/kallisto_linux-v0.50.0.tar.gz -C /opt/biosoft/
echo 'PATH=$PATH:/opt/biosoft/kallisto' >> ~/.bashrc
source ~/.bashrc

# Installing TransDecoder (https://github.com/TransDecoder/TransDecoder/releases)
#wget https://github.com/TransDecoder/TransDecoder/archive/TransDecoder-v5.5.0.tar.gz -P ~/software/
tar zxf ~/software/TransDecoder-v5.5.0.tar.gz -C /opt/biosoft/
mv /opt/biosoft/TransDecoder-TransDecoder-v5.5.0 /opt/biosoft/TransDecoder-v5.5.0
cd /opt/biosoft/TransDecoder-v5.5.0/
make -j 4
echo 'PATH=$PATH:/opt/biosoft/TransDecoder-v5.5.0/' >> ~/.bashrc
source ~/.bashrc

运行数据

# 使用Trinity软件对Illumina数据进行无参考基因组的转录组分析
cd /home/train/08.RNA-seq_analysis_by_trinity
#ln -s ~/03.sequencing_data_preprocessing/?.?.fastq .
# 培训班上笔记本计算性能较差，为了快速计算和减少内存消耗，按如下方法减少数据量再进行计算。
for i in `ls /home/train/03.sequencing_data_preprocessing/?.?.fastq`
do
    x=${i/*\//}
    head -n 1000000 $i > $x
done

# Trinity软件推荐Paired-End Fastq数据的序列名以/1和/2结尾。
# perl -e 'while (<>) { s/^(\@[^\s\/]+).*/$1\/1/; print; $_ = <>; print;  $_ = <>; print; $_ = <>; print; }' read.1.fastq > read_formatted.1.fastq
# perl -e 'while (<>) { s/^(\@[^\s\/]+).*/$1\/2/; print; $_ = <>; print;  $_ = <>; print; $_ = <>; print; }' read.2.fastq > read_formatted.2.fastq

# 使用 Trinity 进行 De novo 组装
Trinity --seqType fq --max_memory 2G --left `ls *.1.fastq | perl -pe 's/\n/,/' | perl -pe 's/,$//'` --right `ls *.2.fastq | perl -pe 's/\n/,/' | perl -pe 's/,$//'` --CPU 8 --jaccard_clip --normalize_reads --output trinity_denovo --bflyCalculateCPU &> trinity_denovo.log
# real  54m34.332s
# user  312m57.809s
# sys   67m20.243s
#由于运行时间过长，可以将其kill掉
[train@MiWiFi-R3P-srv ~]$ kill_procedures_by_name_clf.pl Trinity
# 统计Trinity组装结果
/opt/biosoft/Trinity-v2.15.1/util/TrinityStats.pl trinity_denovo.Trinity.fasta > trinity_denovo.Trinity.fasta.stats
[train@MiWiFi-R3P-srv trinity_denovo]$ cat trinity_denovo.Trinity.fasta.stats 


################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  5566
Total trinity transcripts:  8051
Percent GC: 56.65

########################################
Stats based on ALL transcript contigs:
########################################

    Contig N10: 8102
    Contig N20: 5363
    Contig N30: 4256
    Contig N40: 3424
    Contig N50: 2821

    Median contig length: 1083
    Average contig: 1701.38
    Total assembled bases: 13697830


#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

    Contig N10: 5458
    Contig N20: 3960
    Contig N30: 3135
    Contig N40: 2589
    Contig N50: 2129

    Median contig length: 797.5
    Average contig: 1277.36
    Total assembled bases: 7109784
# 提取最长的Unigene
extract_longest_isoforms_from_TrinityFasta.pl trinity_denovo.Trinity.fasta > trinity_denovo.unigene.longest.fasta

# 使用 Trinity 进行有参考基因组的组装
samtools merge -@ 8 merged.bam ~/06.reads_aligment/hisat2/*.sam
samtools sort -@ 8 -O BAM -o merged.sort.bam merged.bam
Trinity --max_memory 2G --CPU 8 --jaccard_clip --normalize_reads --genome_guided_bam merged.sort.bam --genome_guided_max_intron 4000 --output trinity_genomeGuided --bflyCalculateCPU &> trinity_genomeGuided.log
# real  300m38.653s
# user  676m36.462s
# sys   39m30.345s
/opt/biosoft/Trinity-v2.15.1/util/TrinityStats.pl trinity_genomeGuided/Trinity-GG.fasta > trinity_genomeGuided/Trinity-GG.fasta.stats
#cp ~/08.RNA-seq_analysis_by_trinity/trinity_genomeGuided/Trinity-GG.fasta ~/00.incipient_data/data_for_gene_prediction_and_RNA-seq/
rm merged.*


# 进行表达量计算
mkdir -p /home/train/08.RNA-seq_analysis_by_trinity/exp_cal
cd /home/train/08.RNA-seq_analysis_by_trinity/exp_cal
ln -s ~/00.incipient_data/data_for_gene_prediction_and_RNA-seq/Trinity.fasta ./#转录本序列
/opt/biosoft/Trinity-v2.15.1/util/support_scripts/get_Trinity_gene_to_trans_map.pl Trinity.fasta > Trinity.fasta.gene_trans_map

# 准备转录组序列的索引文件
/opt/biosoft/Trinity-v2.15.1/util/align_and_estimate_abundance.pl --transcripts Trinity.fasta --est_method RSEM --output_dir . --aln_method bowtie2 --prep_reference --gene_trans_map Trinity.fasta.gene_trans_map
# real  0m58.207s
# user  0m57.761s
# sys   0m0.426s

# 进行表达量计算
for i in `ls ../*.1.fastq`
do
    i=${i/*\//}
    i=${i/.1.fastq/}
    echo "/opt/biosoft/Trinity-v2.15.1/util/align_and_estimate_abundance.pl --transcripts Trinity.fasta --seqType fq --left ../$i.1.fastq --right ../$i.2.fastq --est_method RSEM --aln_method bowtie2 --thread_count 8 --gene_trans_map Trinity.fasta.gene_trans_map --output_dir $i &> $i.log"
done > command.align_and_estimate_abundance.list
ParaFly -c command.align_and_estimate_abundance.list -CPU 1 #推荐
#可选sh command.align_and_estimate_abundance.list
# real  32m0.057s
# user  114m33.540s
# sys   8m22.277s
[train@MiWiFi-R3P-srv A]$ lh
total 141M
-rw-r--r-- 1 train train  70M Jul 21 16:48 bowtie2.bam
-rw-r--r-- 1 train train  70M Jul 21 16:48 bowtie2.bam.for_rsem.bam
-rw-r--r-- 1 train train    0 Jul 21 16:48 bowtie2.bam.ok
-rw-r--r-- 1 train train 466K Jul 21 16:49 RSEM.genes.results #基因的表达量
-rw-r--r-- 1 train train 613K Jul 21 16:49 RSEM.isoforms.results #转录本的表达量,TPM可以消除转录本长度不同带来的影响
-rw-r--r-- 1 train train    0 Jul 21 16:49 RSEM.isoforms.results.ok
drwxr-xr-x 2 train train   58 Jul 21 16:49 RSEM.stat
TPM总和100万
补充



for i in `ls */RSEM.isoforms.results`
do
    x=${i/\/RSEM/}
    cp $i $x
done
# cp ~/00.incipient_data/data_for_gene_prediction_and_RNA-seq/RSEM_out/* ./

# 由各个样品的表达量，得到总的表达量矩阵，同时根据表达量去除假阳性序列 
extract_trinity_unigene_by_isofroms_expression_and_ORF.pl Trinity.fasta *.isoforms.results
#计算count总和是否高于100万，如果低于一百万，分析数据不可靠
[train@MiWiFi-R3P-srv exp_cal]$ head out.gene_rawCounts.matrix 
    A   B   C   D   E   F   G
TRINITY_DN0_c0_g1   178.99  190 559.01  71  93  61  68
TRINITY_DN1004_c1_g1    25  22  11  5   5   2   2
TRINITY_DN100_c0_g1 84.99   69  189 20  26  19  13
TRINITY_DN1015_c0_g1    20  25  30  3   0   1   5
TRINITY_DN1035_c0_g1    15  7   39  13  16  17  14
TRINITY_DN1037_c0_g1    50  51  35  15  10  9   9
TRINITY_DN103_c0_g1 56  59  52.01   6   10  6   6
TRINITY_DN1040_c0_g2    40  33  110 19  11  18  20
TRINITY_DN1042_c0_g1    42  49  176 13  16  16  11
[train@MiWiFi-R3P-srv exp_cal]$ cut -f 2 out.gene_rawCounts.matrix | perl -e '<>; while (<>) { $total += $_;} print "$total\n";'
232233.84 #低于100万


# 差异表达分析_样品间的标准化分析,消除样品间的差异
mkdir -p /home/train/08.RNA-seq_analysis_by_trinity/diff_exp
cd /home/train/08.RNA-seq_analysis_by_trinity/diff_exp
ln -s ../exp_cal/out.gene_rawCounts.matrix gene.rawCount.matrix
ln -s ../exp_cal/out.gene_TPM.not_cross.matrix gene.TPM.not_cross.matrix
/opt/biosoft/Trinity-v2.15.1/util/support_scripts/run_TMM_scale_matrix.pl --matrix gene.TPM.not_cross.matrix > gene.TPM.TMM.matrix

echo -e "S4\tA
S4\tB
S4\tC
S2\tD
S2\tE
S2\tF
S2\tG" > samples.txt
# 对两两样本counts总和 >= 40 的基因进行差异表达分析
#/opt/biosoft/trinityrnaseq-Trinity-v2.1.1/Analysis/DifferentialExpression/run_DE_analysis.pl --matrix genes.counts.matrix --method edgeR --samples_file samples.txt --min_rowSum_counts 40 --output edgeR.genes.dir
/opt/biosoft/Trinity-v2.15.1/Analysis/DifferentialExpression/run_DE_analysis.pl --matrix gene.rawCount.matrix --method edgeR --samples_file samples.txt --output edgeR.genes.dir

cd edgeR.genes.dir
[train@MiWiFi-R3P-srv edgeR.genes.dir]$ ls
gene.rawCount.matrix.S2_vs_S4.edgeR.count_matrix  gene.rawCount.matrix.S2_vs_S4.edgeR.DE_results.MA_n_Volcano.pdf
gene.rawCount.matrix.S2_vs_S4.edgeR.DE_results    gene.rawCount.matrix.S2_vs_S4.S2.vs.S4.EdgeR.Rscript
[train@MiWiFi-R3P-srv edgeR.genes.dir]$ head gene.rawCount.matrix.S2_vs_S4.edgeR.DE_results
sampleA sampleB logFC   logCPM  PValue  FDR
TRINITY_DN190_c0_g1 S2  S4  8.99472274370147    15.8184438491505    2.64453178222291e-123   1.81943786616937e-120
TRINITY_DN1984_c0_g1    S2  S4  7.77684165600643    13.7102324085575    1.8328163121404e-117    6.30488811376297e-115
TRINITY_DN14_c0_g1  S2  S4  7.01104197690975    18.0524459121504    6.43603234776519e-93    1.47599675175415e-90
TRINITY_DN1460_c0_g1    S2  S4  6.55399916827941    8.70473121995052    1.3399809141212e-65 2.30476717228846e-63
TRINITY_DN198_c0_g1 S2  S4  3.87279572531066    10.0995425444949    1.04447111920147e-59    1.43719226002122e-57
TRINITY_DN400_c0_g1 S2  S4  7.18963166670787    8.10209691711092    1.34154919588681e-57    1.53830974461688e-55
TRINITY_DN4100_c0_g1    S2  S4  7.49481831648387    8.10572476735801    6.0425927424905e-57 5.93900543833352e-55
TRINITY_DN96_c0_g1  S2  S4  5.82899775763059    12.2441502545614    4.70019463324168e-54    4.04216738458785e-52
TRINITY_DN221_c0_g1 S2  S4  3.83239810737771    8.84603355603372    1.18232923097752e-52    9.03825012125039e-51
#一般pvalue和FDR，参考FDR
# 提取差异表达基因进行聚类分析和热图制作
/opt/biosoft/Trinity-v2.15.1/Analysis/DifferentialExpression/analyze_diff_expr.pl --matrix ../gene.TPM.TMM.matrix -P 0.001 -C 2 --samples ../samples.txt
#/opt/biosoft/Trinity-v2.15.1/Analysis/DifferentialExpression/analyze_diff_expr.pl --matrix ../gene.TPM.TMM.matrix -P 0.01 -C 1 --samples ../samples.txt --order_columns_by_samples_file

# 根据距离结果进行分类
# 自动分类
# 共表达趋势分析
/opt/biosoft/Trinity-v2.15.1/Analysis/DifferentialExpression/define_clusters_by_cutting_tree.pl --Ptree 50 -R diffExpr.P0.001_C2.matrix.RData
# 手动分类
# R
# > load("diffExpr.P0.001_C2.matrix.RData")
# > source("/opt/biosoft/Trinity-v2.15.1/Analysis/DifferentialExpression/R/manually_define_clusters.R")
# > manually_define_clusters(hc_genes, data)
# 然后左键点击选择子类，右键结束选择。
# 结果生成了文件夹 manually_defined_clusters_4 。其中 4 表示手动分成了 4 类。
# cd manually_defined_clusters_4
# 对每类结果进行趋势线的绘制
# /opt/biosoft/Trinity-v2.15.1/Analysis/DifferentialExpression/plot_expression_patterns.pl cluster_*

image.png


sudo cpan -i Statistics::Basic
DEG_compare_2_samples.pl
Usage:
    /home/train/bin/DEG_compare_2_samples.pl [options] --FDR_edgeR 0.001 --FDR_DESeq2 0.001 --log2FC 2 genes.TPM_TMM.matrix Sample1:LabelA,LabelB,LabelC Sample2:LabelD,LabelE,LabelF

    程序用于分析两个样本的差异表达基因。注意事项：
    （1）程序输入的第一个信息为所有样本的表达量矩阵文件，第二个信息为样本1的标签信息，第二个信息为样本2的标签信息。样本和标签之间用冒号分割，标签之间用逗号分割。
    （2）输入的表达量矩阵文件中第一行一定要包含后面输入的所有Label名称，Label名称不推荐含有怪异字符，特别是逗号和冒号等。表达量矩阵文件中可以包含有更多Label的表达量信息，但程序会忽略之。
    （3）表达量矩阵文件中的数据需要为使用TMM或其它算法进行了样本间标准化的表达量。程序进行edgeR分析时不会额外再进行TMM标准化了。推荐使用normalizing_TPM_matrix_by_TMM.pl命令对TPM表达量进行TMM标准化。
    （4）程序仅读取后面的两个样本信息，样本名字可以随意取，但推荐不要含有怪异字符。
    （5）--FDR_edgeR、--FDR_DESeq2和--log2FC参数必须至少有一个，或推荐全部都有。
    （6）程序输出三个结果文件，DE_result.tab包含所有基因的FDR、log2FC和表达量信息，Up.tab和Down.tab分别是在第二个样本中上调或下调的基因信息，当使用--annot_file参数提供注释文件时，在第二列插入基因的功能注释信息。三个输出文件均按FDR和log2FC进行了排序。
    （7）程序在标准错误输出中给出用于差异分析的最小表达量值、各个算法得到的差异表达基因数量。

    参数：
    --tmp_dir <string>    default: tmp_$date_$pid
    程序运行时临时文件夹名称。

    --out_prefix <string>    default: Sample1_VS_Sample2
    设置输出文件前缀。

    --FDR_edgeR <float>    default: None
    设置该参数后，则使用edgeR进行差异表达基因分析，并要求差异表达量基因的FDR <= 该阈值，推荐设置为 0.001。若不设置该参数值，则程序不会调用edgeR进行分析。

    --edgeR_dispersion <float>    default: 0.1
    若两个样本都没有生物学重复，使用edgeR的exactTest进行p值计算时，设置其dispersion参数的值。

    --FDR_DESeq2 <float>    default: None
    设置该参数后，则使用DESeq2进行差异表达基因分析，并要求差异表达量基因的FDR <= 该阈值，推荐设置为 0.001。若不设置该参数值，则程序不会调用DEseq2进行分析。当两个样本都没有生物学重复时，不会进行DESeq2分析。虽然无生物学重复DESeq2运行不报错，但是一般得不到差异基因，导致分析没有意义。

    --log2FC <float>    default: None
    设置该参数后，则使用表达量差异倍数进行差异基因分析，并要求差异表达量基因的 |log2表达量差异倍数| >= 该阈值，推荐设置为 2。若不设置该参数值，则程序不会则使用表达量差异倍数进行分析。

    --annot_file <string>    default: None
    若输入基因功能注释文件，则会在输出文件中包含基因的功能注释信息。要求文件使用制表符分割，第一列是基因ID，后面的一列或多列注释信息都能被识别，也支持一个基因在多行上有注释。程序会将一个基因所有行所有列的注释信息用 | 符号分割，赋值为其基因的注释信息，放入到差异表达基因结果文件的第二列。

    --min_value <float>    default: None
    设置最小的表达量值。当基因的表达量低于该值时，则将其表达量修改为该值，再用于差异基因分析。

    --median_ratio <float>    default: 0.01
    若未设置--min_value参数值，则使用该参数值按如下步骤计算出最小表达量值：（1）去除表达量矩阵中所有不为0的数据，求中位数；（2）再去除表达量值 < 中位数 * 该参数值的数据，再次求中位数；（3）使用第二次的中位数 * 该参数设定值即得到最小表达量值。

    --keep_tmp    default: None
    添加该参数，程序会保留临时文件夹。

    --help    default: None
    display this help and exit.

功能注释

[train@MiWiFi-R3P-srv transdecoder]$ ls ../exp_cal/out.*
../exp_cal/out.effective.fasta  #转录本序列文件，文件较大
../exp_cal/out.gene_rawCounts.matrix  
../exp_cal/out.gene_TPM.not_cross.matrix  
../exp_cal/out.unigene.fasta #unigene文件，使用该文件

# 使用 transdecoder 进行ORF预测
mkdir -p /home/train/08.RNA-seq_analysis_by_trinity/transdecoder
cd /home/train/08.RNA-seq_analysis_by_trinity/transdecoder
#ln -s ../exp_cal/out.unigene.fasta Unigene.fasta
cp ~/00.incipient_data/data_for_gene_prediction_and_RNA-seq/Unigene.fasta ./

# 首先，根据6个读码框翻译得到蛋白序列
TransDecoder.LongOrfs -t Unigene.fasta -m 100 #-m 100 至少有100个氨基酸
# real  0m11.108s
# user  0m11.068s
# sys   0m0.042s
perl -pe 's/^>(\S+).*/>$1/' Unigene.fasta.transdecoder_dir/longest_orfs.pep > longest_orfs.pep #寻找每个转录本最长的开放性阅读框

# 然后，将ORF蛋白序列比对到公共数据库，能比对上的则表示存在对应转录本序列
# 比对到SwissProt数据库
#wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz -P ~/software/
#gzip -dc ~/software/uniprot_sprot.fasta.gz > uniprot_sprot.fasta
#diamond makedb --threads 8 --db uniprot_sprot --in uniprot_sprot.fasta
#diamond blastp --db uniprot_sprot --query longest_orfs.pep --out blast.xml --outfmt 5 --sensitive --max-target-seqs 20 --evalue 1e-5 --id 20 --tmpdir /dev/shm --index-chunks 1
#parsing_blast_result.pl --no-header --max-hit-num 20 --query-coverage 0.2 --subject-coverage 0.1 blast.xml > blastp.outfmt6
cp ~/00.incipient_data/data_for_gene_prediction_and_RNA-seq/blastp.outfmt6 ./

# 比对到PFAM数据库
#para_hmmscan.pl --cpu 40 --chunk 30 --no_cut_ga --hmm_db /opt/biosoft/bioinfomatics_databases/Pfam/Pfam-AB.hmm longest_orfs.pep | grep -P "^#" -v > pfam.domtbl
cp ~/00.incipient_data/data_for_gene_prediction_and_RNA-seq/pfam.domtbl ./

# 最后，去除假阳性ORF，保留正确的ORF
TransDecoder.Predict -t Unigene.fasta --retain_pfam_hits pfam.domtbl --retain_blastp_hits blastp.outfmt6
# real  0m51.079s
# user  0m50.418s
# sys   0m0.671s
#统计ORF名称
[train@MiWiFi-R3P-srv transdecoder]$ grep ">" Unigene.fasta.transdecoder.pep  | perl -e 'while (<>) { print "$1\n" if m/^>(\w+)/;}' | head
TRINITY_DN0_c0_g1
TRINITY_DN0_c0_g1
TRINITY_DN0_c0_g1
TRINITY_DN0_c0_g1
TRINITY_DN0_c1_g1
TRINITY_DN0_c1_g1
TRINITY_DN1000_c0_g2
TRINITY_DN1000_c0_g2
TRINITY_DN1002_c0_g1
TRINITY_DN1004_c0_g1
[train@MiWiFi-R3P-srv transdecoder]$ grep ">" Unigene.fasta.transdecoder.pep  | perl -e 'while (<>) { print "$1\n" if m/^>(\w+)/;}' | uniq | wc -l
2266




# 检查正义链负义链的ORF长度和数量
perl -e 'while (<>) { if (m/\tCDS\t(\d+)\t(\d+)\t\S+\t(\S+)\t/) { $length = $2 - $1 + 1; push @stats1, $length if $3 eq "+"; push @stats2, $length if $3 eq "-"; } } @stats1 = sort {$a <=> $b} @stats1; @stats2 = sort {$a <=> $b} @stats2; $stats1 = @stats1; $stats2 = @stats2; print "strand\tnum\tmedian_length\n+\t$stats1\t$stats1[$stats1/2]\n-\t$stats2\t$stats2[$stats2/2]\n";' Unigene.fasta.transdecoder.gff3

# 测试连特异性测序转录组De novo组装序列
#Trinity --seqType fq --max_memory 2G --left ../E.1.fastq --right ../E.2.fastq --CPU 40 --jaccard_clip --normalize_reads --output trinity_denovo_SS --bflyCalculateCPU --SS_lib_type RF &> trinity_denovo_SS.log
#extract_longest_isoforms_from_TrinityFasta.pl trinity_denovo_SS.Trinity.fasta > trinity_denovo_SS.unigene.longest.fasta
cp ~/00.incipient_data/data_for_gene_prediction_and_RNA-seq/trinity_denovo_SS.unigene.longest.fasta ./
TransDecoder.LongOrfs -t trinity_denovo_SS.unigene.longest.fasta -m 100
TransDecoder.Predict -t trinity_denovo_SS.unigene.longest.fasta
# 检测在正负链上同时进行预测ORF情况
perl -e 'while (<>) { if (m/\tCDS\t(\d+)\t(\d+)\t\S+\t(\S+)\t/) { $length = $2 - $1 + 1; push @stats1, $length if $3 eq "+"; push @stats2, $length if $3 eq "-"; } } @stats1 = sort {$a <=> $b} @stats1; @stats2 = sort {$a <=> $b} @stats2; $stats1 = @stats1; $stats2 = @stats2; print "strand\tnum\tmedian_length\n+\t$stats1\t$stats1[$stats1/2]\n-\t$stats2\t$stats2[$stats2/2]\n";' trinity_denovo_SS.unigene.longest.fasta.transdecoder.gff3
# 检测仅在正以链上进行预测ORF情况
perl -e 'while (<>) { if (m/\tCDS\t(\d+)\t(\d+)\t\S+\t(\S+)\t/) { $length = $2 - $1 + 1; push @stats1, $length if $3 eq "+"; push @stats2, $length if $3 eq "-"; } } @stats1 = sort {$a <=> $b} @stats1; @stats2 = sort {$a <=> $b} @stats2; $stats1 = @stats1; $stats2 = @stats2; print "strand\tnum\tmedian_length\n+\t$stats1\t$stats1[$stats1/2]\n-\t$stats2\t$stats2[$stats2/2]\n";' ~/00.incipient_data/data_for_gene_prediction_and_RNA-seq/trinity_denovo_SS.unigene.longest.fasta.transdecoder.gff3
cd ..

组装后的文件就是转录本文件
在未消除样品之间的差异之前，用TPM来比较勉强可以，但是FPKM不行。

转录组入门学习（四）
今天进行序列比对课程的学习序列比对 1. 无参转录组 1.1 使用拼接工具组装转录本trinity 1.2 tr...
使用Trinity 软件进行拼接（无参转录组）
文章背景（此背景非彼背景）本文主要是对Trinity软件的背景、安装以及参数做一个详细的介绍。别问为什么没有...
全基因组测序从头测序(de novo sequencing)
全基因组测序全基因组测序分为从头测序(de novo sequencing)和重测序(re-sequencing...
转录组数据库的基本使用（一）-GO数据库
转录组数据库的基本使用（一）-GO数据库转录组数据库有哪些？目前在对转录组数据进行分析的时候，很多测序公司通常...
转录组组装软件：Trinity (转载)
Trinity是由 the Broad Institute 开发的转录组denovo组装软件，由三个独立的软件模块...
空间转录组分析-工具箱
空间转录组注释通过marker进行spot注释 1、手动2、使用工具-SingleR空间转录组第九讲：利用Sin...
转录组测序技术和结果解读（五）——参考基因组mapping
转录组测序根据目标物种是否存在基因组信息分为有参转录组和无参转录组。有参转录组通过将测序数据与参考基因组进行ma...
基因组从头组装软件
从头组装- de novo assembly 指在不需要任何参考序列的情况下对某一物种进行基因组测序，然后将测得的...
基因组结构注释软件列表
基因组结构注释主要包括三个方面：从头注释（de novo prediction）、同源注释（homology-ba...
4、RNAseq(4)--Hisat2进行序列比对及Samtoo
概况：使用处理后的fastq文件和基因组与转录组比对，确定在转录组或者基因组中的关系。在转录组和基因组的比对采取的...