美文网首页组学学习
跟着Nature Communications学习Hisat-T

跟着Nature Communications学习Hisat-T

作者: 小杜的生信筆記 | 来源:发表于2023-12-02 18:05 被阅读0次

    一边学习,一边总结,一边分享!

    详细教程请访问:
    组学分析流程

    本期分析流程

    1. Hisat2-Samtools
    2. Trinity_GG_denovo
    3. 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的分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!!

    相关文章

      网友评论

        本文标题:跟着Nature Communications学习Hisat-T

        本文链接:https://www.haomeiwen.com/subject/lnghgdtx.html