GATK参数报错求助：批量BAM文件局部重排代码修正

Hey fellow Bash/GATK developer! Let's work through the issues in your script and get that local realignment running smoothly.

Key Problems in Your Original Code #

• You left the -I (input BAM) parameter incomplete — it needs to point to the file you're looping over
• Using ls to iterate over files can break if filenames have spaces or special characters (we'll swap this for a safer glob pattern)
• GATK3's IndelRealigner relies on two critical components you're missing:
  1. A target intervals file (generated first via RealignerTargetCreator)
  2. Known indel reference data (passed with the -known flag)

Corrected Complete Workflow #

Local realignment with GATK3 is a two-step process: first identify regions that need realignment, then perform the realignment itself.

Step 1: Generate Realignment Target Intervals

# First, create target intervals for each BAM file
echo "Generating realignment target intervals..."
for file in adp/map/*marked_duplicates.bam; do
  # Create a consistent output filename for intervals
  target_file="${file%.bam}.intervals"
  
  java -jar ~/software/GenomeAnalysisTK-3.3-0/GenomeAnalysisTK.jar \
    -T RealignerTargetCreator \
    -R ~/flybase/fb-r5.57/dmel-all-chromosome-r5.57.fasta \
    -I "$file" \
    -known /path/to/your/drosophila_known_indels.vcf \  # Replace with your actual known indels file path
    -o "$target_file"
  
  echo "Created target intervals: $target_file"
done

Step 2: Run Indel Realignment

# Now execute the local realignment using the generated targets
echo "Starting local realignment..."
for file in adp/map/*marked_duplicates.bam; do
  # Define output filename for realigned BAM
  realigned_bam="${file%.marked_duplicates.bam}.realigned.bam"
  target_file="${file%.bam}.intervals"
  
  java -jar ~/software/GenomeAnalysisTK-3.3-0/GenomeAnalysisTK.jar \
    -T IndelRealigner \
    -R ~/flybase/fb-r5.57/dmel-all-chromosome-r5.57.fasta \
    -I "$file" \
    -targetIntervals "$target_file" \
    -known /path/to/your/drosophila_known_indels.vcf \  # Same known indels file as above
    -o "$realigned_bam"
  
  echo "Completed realignment for: $realigned_bam"
done

Quick Tips #

• Replace /path/to/your/drosophila_known_indels.vcf with the actual path to your FlyBase indel reference file
• Wrapping variables like $file in double quotes prevents bugs with spaces or special characters in filenames
• The parameter expansion ${file%.bam} lets us generate clean, consistent output filenames without repeating the original path manually

内容的提问来源于stack exchange，提问作者multiverse