#run hisat2 for gall vs leaf RNAseq #trim reads mkdir trimmed_reads cd trimmed_reads/ for i in ../rawdata/*_1.fastq.gz do file=`basename $i` numbertwo=${i//_1.fastq.gz/_2.fastq.gz} numbertwofile=${file//_1.fastq.gz/_2.fastq.gz} bsub -n 2 -o gall_trim_reads -We 60 "cutadapt -j 2 --trim-n -a revcompTruSeq_Univ_Adap=AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -A revcompTruSeq_Univ_Adap=AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \ -a TruSeq_Idx_Adap=GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG -A TruSeq_Idx_Adap=GATCGGAAGAGCACACGTCTGAACTCCAGTCACNNNNNNATCTCGTATGCCGTCTTCTGCTTG -m 36 -q 10,10 --nextseq-trim=10 -o catrim.gall.${file} -p \ catrim.gall.${numbertwofile} ${i} ${numbertwo}" done conda activate RNAseq ln -s ~/Hvir_referenceGenome/Hvir_nuclei_sn_run2_2_pseudohap_mskd.smplchr.fa . ln -s ~/Hvir_referenceGenome/UTRme/Hvir-UTRme_v1-3-best-score_cleanUTRs_v2.gff . #bsub -n 8 -o hisat2_idx "hisat2-build ~/Hvir_referenceGenome/Hvir_nuclei_sn_run2_2_pseudohap_mskd.smplchr.fa \ # Hvir_nuclei_sn_run2_2_pseudohap_mskd.smplchr_hisat2 -p 8" #map with hisat2 cd trimmed_reads for i in catrim*1.fastq.gz do root_name=`basename ${i} _1.fastq.gz` second_name=${root_name}_2.fastq.gz bsub -n 4 -o hisat2_aln_${root_name} \ "hisat2 -x ~/Hvir_RNAseq/RedVsGreen_Sep2020_UTR_hisat2/Hvir_nuclei_sn_run2_2_pseudohap_mskd.smplchr_hisat2 \ -p 4 \ -1 ${i} -2 ${second_name} | samtools sort -O BAM -@ 4 - > ${root_name}.bam" done cd .. #reformat gff grep exon Hvir-UTRme_v1-3-best-score_cleanUTRs_v2.gff | sed 's|Parent|transcript_id|' > exons.gff3 sed 's|.*transcript_id=\(g.*\).t.*|; gene_id=\1|' exons.gff3 > genes.gff3 paste -d \0 exons.gff3 genes.gff3 > Hvir-UTRme_v1-3-best-score_cleanUTRs_v2.exons.gff3 #htseq to generate counts mkdir counts for i in trimmed_reads/*.bam do root_name=`basename ${i} .bam` prefix=${root_name//catrim./} bsub -o counts/count.${prefix}.txt \ "htseq-count -i "gene_id" --stranded=no -f bam -t exon \ ${i} Hvir-UTRme_v1-3-best-score_cleanUTRs_v2.exons.gff3 > counts/${prefix}.counts" done conda deactivate