跟着Nature Communications学习Hisat-Trinity-PASA等分析流程
2023-12-14 17:41:30
一边学习,一边总结,一边分享!
详细教程请访问:
组学分析流程
本期分析流程
- Hisat2-Samtools
- Trinity_GG_denovo
- PASA
…
本期教程文章
题目:Genomic insights into local adaptation and future climate-induced vulnerability of a keystone forest tree in East Asia
Hisat2-samtools分析流程
#!/bin/bash
genome=$1
index=${genome%.*}
rna_1_fq=`cat $2|grep 1P|sed ":a;N;s/\n/,/g;ta"` #1.fq path list
rna_2_fq=`cat $2|grep 2P|sed ":a;N;s/\n/,/g;ta"` #2.fq path list
#echo $index
hisat2-build -p 20 $genome $index
hisat2 -x $index \
-1 $rna_1_fq\
-2 $rna_2_fq\
--threads 20 \
--min-intronlen 20 \
--max-intronlen 20000 \
--dta \
--score-min L,0.0,-0.4 \
-S ${index}.sam
samtools sort -@ 20 \
-o ${index}.sorted.bam \
-O BAM \
${index}.sam
PSSA_align
#!/bin/bash
export PATH="$PATH:/usr_storage/jcf/.conda/envs/PASA"
source /pub_storage2/new_PASA/.bashrc
#cat $Trinity_GG $Trinity_denovo >transcripts.fasta #
transcripts_fasta="$1" # transcripts.fasta generated from merging fasta file of Trinity denovo and Trinity genome guided mode
#perl -e 'while(<>) { print "$1\n" if />(\S+)/ }' Trinity.fasta >tdn.accs #
denovo_transcript_id="$2"
alignAssembly_config="$3"
genome="$4" #reference fasta file
seqclean $transcripts_fasta \
-v /pub_storage2/PASA/UniVec
Launch_PASA_pipeline.pl -c $alignAssembly_config \
-C -R -T \
-g $genome \
-t $transcripts_fasta.clean \
-u ${transcripts_fasta} \
--ALIGNERS gmap,blat \
--CPU 8 \
--TDN $denovo_transcript_id
Trinity GG denovo
#!/bin/bash
#conda activate trinity
export PATH="$PATH:/usr_storage/jcf/.conda/envs/trinity"
rna_1_fq="cat $1|sed ":a;N;s/\n/,/g;ta"" #1.fq path list
rna_2_fq="cat $2|sed ":a;N;s/\n/,/g;ta"" #2.fq path list
bam="$3" #sorted.bam from hisat
out=${bam%.*}
Trinity --left $rna_1_fq \
--right $rna_2_fq \
--seqType fq \
--max_memory 100G \
--no_normalize_reads \
--CPU 20 \
--bflyCalculateCPU \
--output trinity_denovo_$out
Trinity --genome_guided_bam $bam \
--genome_guided_max_intron 10000 \
--max_memory 100G \
--no_normalize_reads \
--CPU 20 \
--bflyCalculateCPU\
--output trinity_GG_$out
ab homo
#!/bin/bash
export PATH="$PATH:/usr_storage/jcf/.conda/envs/BUSCO"
source /usr_storage/jcf/geta-user204/.bashrc
rna_1_fq="cat $1|sed ":a;N;s/\n/,/g;ta"" #1.fq path list
rna_2_fq="cat $2|sed ":a;N;s/\n/,/g;ta"" #2.fq path list
genome="$3" #genome fasta file
conf="$4" #small genome conf.txt of geta pipepline setting as default parameters
out=${genome%.*}
homo_pro="$5"
geta.pl \
--RM_species Embryophyta\
--out_prefix `pwd`/$out \
--config $conf \
--cpu 20 \
--protein $homo_pro\
-genome $genome \
-1 $rna_1_fq \
-2 $rna_2_fq \
--augustus_species $out
Evm
#!/bin/bash
export PATH="/usr_storage/xyf/jcf/genewise/EVM/EVidenceModeler-1.1.1/EvmUtils/:$PATH"
genome="$1" #genome fasta file
augustus_gff3="$2" #gff3 generated from augutus
genewise_gff3="$3" #gff3 generated from tblastn and genewise
pasa_align_gff3="$4" #gff3 generated from PASA
repeat_gff3="$5" #repeat gff3 generated from repeatemasker
partition="$6" #partition path for evm
partition_EVM_inputs.pl \
--genome $genome\
--gene_predictions $augustus_gff3 \
--protein_alignments $genewise_gff3 \
--transcript_alignments $pasa_align_gff3 \
--repeats $repeat_gff3 \
--segmentSize 5000000 \
--overlapSize 10000 \
--partition_listing $partition
write_EVM_commands.pl \
--genome $genome \
--gene_predictions $augustus_gff3 \
--protein_alignments $genewise_gff3 \
--transcript_alignments $pasa_align_gff3 \
--repeats $repeat_gff3 \
--output_file_name evm.out \
--weights $weight >command.list
ParaFly -c command.list -CPU 32
recombine_EVM_partial_outputs.pl \
--partitions $partition \
--output_file_name evm.out
convert_EVM_outputs_to_GFF3.pl \
--partitions $partition \
--output_file_name evm.out \
--genome $genome
cat */evm.out.gff3 >evm.out.gff3
PASA update
#!/bin/bash
export PATH="$PATH:/usr_storage/jcf/.conda/envs/PASA "
source /pub_storage2/new_PASA/.bashrc
genome="$1" #genome fasta file
annotation_conf="$2" #pasa annotation compare conf
transcripts_fasta="$3" #transcripts_fasta file for PASA seqclean step
gff3="$4" #gff3 for PASA updata
Launch_PASA_pipeline.pl \
-c $annotation_conf\
-A -T -L \
-g $genome\
-t ${transcripts_fasta}.clean \
-u $transcripts_fasta \
--annots $gff3
这里只是提供了各个分析流程的脚本,对于初学者来说是比较有好的。我们在转录组上游分析教程[零基础]中提供了详细转录组上游分析的参数,对于初学者来说是比较友好的。
往期文章:
1. 复现SCI文章系列专栏
2. 《生信知识库订阅须知》,同步更新,易于搜索与管理。
3. 最全WGCNA教程(替换数据即可出全部结果与图形)
4. 精美图形绘制教程
5. 转录组分析教程
小杜的生信筆記 ,主要发表或收录生物信息学的教程,以及基于R的分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!!
文章来源:https://blog.csdn.net/kanghua_du/article/details/134767951
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符,请联系我的编程经验分享网邮箱:veading@qq.com进行投诉反馈,一经查实,立即删除!
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符,请联系我的编程经验分享网邮箱:veading@qq.com进行投诉反馈,一经查实,立即删除!