https://bioinformaticsonline.com/snippets/view/44725/bash-script-to-discover-pirna-in-transcriptome-data? https://bioinformaticsonline.com/snippets/view/44725/bash-script-to-discover-pirna-in-transcriptome-data Fri, 13 Dec 2024 11:47:00 -0600 https://bioinformaticsonline.com/snippets/view/44725/bash-script-to-discover-pirna-in-transcriptome-data #!/bin/bash # Variables (modify these as per your setup) INPUT_FASTQ="input_reads.fastq" ADAPTER_SEQ="TGGAATTCTCGGGTGCCAAGG" REFERENCE_GENOME="reference_genome.fa" BOWTIE_INDEX="reference_index" OUTPUT_DIR="piRNA_analysis" THREADS=4 # Create output directory mkdir -p $OUTPUT_DIR # Step 1: Quality Control echo "Running FastQC for quality control..." fastqc $INPUT_FASTQ -o $OUTPUT_DIR # Step 2: Adapter Trimming echo "Trimming adapters with Cutadapt..." cutadapt -a $ADAPTER_SEQ -o $OUTPUT_DIR/trimmed_reads.fastq $INPUT_FASTQ # Step 3: Mapping Reads to Reference Genome echo "Mapping reads to reference genome using Bowtie..." bowtie -v 1 -k 1 --best -p $THREADS $BOWTIE_INDEX $OUTPUT_DIR/trimmed_reads.fastq -S $OUTPUT_DIR/aligned_reads.sam # Step 4: Convert SAM to BAM and Sort echo "Converting SAM to BAM and sorting..." samtools view -Sb $OUTPUT_DIR/aligned_reads.sam | samtools sort -o $OUTPUT_DIR/sorted_reads.bam # Step 5: Extract Reads of piRNA Size (24-32 nt) echo "Filtering reads of size 24-32 nt..." bedtools bamtofastq -i $OUTPUT_DIR/sorted_reads.bam -fq $OUTPUT_DIR/all_reads.fastq seqkit seq -m 24 -M 32 $OUTPUT_DIR/all_reads.fastq > $OUTPUT_DIR/piRNA_size_reads.fastq # Step 6: Detect Sequence Bias (Optional) echo "Checking sequence bias using WebLogo-compatible data..." seqkit fx2tab $OUTPUT_DIR/piRNA_size_reads.fastq | cut -f2 | awk '{print ">seq"NR"\n"$0}' > $OUTPUT_DIR/piRNA_sequences.fa # Step 7: Identify piRNA Clusters # This step requires a tool like proTRAC or PIRANHA. Example placeholder: echo "Identifying piRNA clusters (requires proTRAC or PIRANHA)..." # Example with proTRAC: # proTRAC.pl -s $OUTPUT_DIR/sorted_reads.bam -g $REFERENCE_GENOME -o $OUTPUT_DIR/clusters # Step 8: Annotate Clusters # Annotation depends on your genome's annotation file # bedtools intersect example placeholder: # bedtools intersect -a clusters.bed -b genome_annotation.gtf > annotated_clusters.bed # Step 9: Clean up intermediate files (optional) echo "Cleaning up intermediate files..." rm $OUTPUT_DIR/aligned_reads.sam $OUTPUT_DIR/all_reads.fastq # Done echo "piRNA discovery pipeline completed! Results are in $OUTPUT_DIR."]]> Abhi