ENCODE dynamic k-mer–assisted mapping

1. Paired-end reads are mapped using bowtie2 with the following parameters:
   1. --no-discordant
   2. --no-mixed
   3. --very-sensitive
   4. --no-unal
   5. --omit-sec-seq
   6. --xeq
   7. --reorder
2. Reads are filtered using samtools view and the parameters "-b -F 1804 -f 2 -q 2". This removes reads with a quality score of 1 (multi-mapping, as marked by bowtie2), reads that do not have the "properly aligned" flag (bit 2), and reads that have any of the following flags: unmapped segment (bit 4), unmapped next segment (bit 8), secondary alignment (bit 256), failed filter (bit 512), and duplicate (bit 1024).
3. Duplicate reads are removed using picard markduplicates and the parameters "VALIDATION_STRINGENCY=LENIENT ASSUME_SORT_ORDER=queryname REMOVE_DUPLICATES=true". This removes reads originating from the same molecule as determined by comparison of their 5' positions and sequences. The read with the highest base quality sum is retained from the pool of duplicate reads.
4. Reads from replicate samples are pooled using samtools merge.
5. Pooled reads are sorted by name. This places read ends next to each other in the order.
6. K-mer databases are created for the genome using KMC2 and parameters "-k<KMERSIZE> -fm -ci2" with KMER_SIZE taking the value 50, 75, 80, 85, 90, 95, and 100 for a total of 7 k-mer databases, The remaining parameters specify that there are multiple sequences in the genome fasta file and to exclude k-mers that only appear once (unique) in the genome.
7. For each k-mer database, a count map wiggle file is produced using kmc_genome_counts and the genome fasta file. This wiggle file has a value for each base position specifying the number of times the corresponding k-mer from the genomic sequence appears anywhere in the genome. The k-mer sequence starts at the base position and extend downstream a number of bases equal to the k-mer size.
8. The k-mer count map wiggle files are converted to bigwig files using the UCSC genome browser utility wigToBigWig.
9. Pooled and sorted reads are filtered by unique k-mers using the custom script filter_by_unique_kmers.py found at https://github.com/msauria/T2T_Encode_Analysis/blob/main/bin/filter_by_unique_kmers.py and the k-mer bigwig files. This script loads in each of the k-mer count map bigwigs (specified by a filename template and comma-separated list of k-mer sizes). For each read, the reference stop position minus the reference start position (the reference span) is used to determine the k-mer length to use, selecting the largest k-mer size that is equal to or less than the reference span. Once the k-mer size has been determined, the score for each k-mer of that size occurring in the reference sequence covered by the mapped read is checked in the corresponding k-mer count map (e.g. a read mapping to 99 bases of reference sequence would use the k-mer size 95, and would contain 5 separate 95-mers that would be looked up in the count map). If any of the values in the count map corresponding to the contained k-mers is equal to one (signifying a unique k-mer sequence), the read and its mate are retained. Otherwise the read and its mate are discarded.
10. The k-mer-filtered reads are sorted by position using samtools sort.

Software versions and availability:

• Bowtie2 - v0.69.0
• Samtools - v1.10
• Picard - v0.68.0
• KMC - https://github.com/msauria/KMC/tree/kmer_mapping extended from v3.1.2rc1
• kmc_genome_count - included in the above KMC package
• wigToBigWig - v4
• filter_by_unique_kmers.py - https://github.com/msauria/T2T_Encode_Analysis, requires:
  • pysam - v0.16.0
  • pyBigWig - part of deeptools v3.4.3
  • numpy - v1.19.2