Indexing fragments file created using a custom script (paired RNA, ATAC single-end reads)

[sallyseullee-0821]

Hi,

I used custom scripts to process paired RNA and ATAC data (single-end reads) into matrix files and bam output files for joint analysis, and also fragments.tsv.gz for the ATAC data.
I am trying to create ATAC assay using CreateChromatinAssay with the fragments file I generated. However, it seems like I require an indexed fragments file in the same directory.

The code I used to generate my fragments.tsv.gz from the bam files is the following:

my $bam = $ARGV[0]; #or die $!;
#### usage perl bam2fragments.pl input.bam

open OUT, "|gzip - > $ARGV[0]\_fragments.tsv.gz" or die $!;

# chr1	3000050	3000595 Barcode	1

open IN, "samtools view $bam|" or die $!;
while(<IN>){
  chomp;
  my @tmp = split/\s+/, $_;
  my $bc = substr($tmp[0], -19, 8);
  my $chr = $tmp[2];
  my $pss = $tmp[3];
  my $pse = $pss + 150;
  print OUT "$chr\t$pss\t$pse\t$bc\t1\n";

}
close IN;
close OUT;

With the generated tsv.gz file, I used

gzip -d <fragments>
bgzip <fragments>
tabix -p bed <fragments>

However, I get an error saying

[E::hts_idx_push] Chromosome blocks not continuous.

I would really appreciate your help with troubleshooting this.
Thank you! Looking forward to any insights!

Best,
Sally

[timoast]

The error indicates that the fragments are not correctly sorted by coordinate.

FYI, the fragment file you're constructing here will not be correct. The start and end coordinates of each entry represent the two Tn5 cut sites that generated the fragment. In this case, you have one cut site and then are adding another one 150 bp away which is not necessarily where Tn5 cut.

[sallyseullee-0821]

Thank you for your prompt comment! And you are right; since the paired datasets I am working with are single-end sequenced, I do not have the fragment length information. So, I used 151 bp for all of them.

Then, does this workflow seem correct?

gzip -d fragments.tsv.gz
sort -k1,1 -k2,2n -k3,3n fragments.tsv
bgzip fragments.tsv
tabix -p bed fragments.tsv.gz

[timoast]

I don't think you need -k3:

sort -k 1,1 -k2,2n frags.bed > frags.sort.bed
bgzip frags.sort.bed
tabix -p bed frags.sort.bed

[timoast]

Also, the PCR duplicate collapsing is not being done here...each read would be considered a different fragment. Take a look here:
https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/output/fragments

Your results using this file would be inaccurate due to incorrect placing of the second Tn5 cut site and lack of PCR duplicate removal.