Pysam

Lecture d'un fichier (ne nécessite pas d'index) :

reader = pysam.AlignmentFile('myFile.bam', 'rb')
for align in reader.fetch(until_eof = True):
  print(align.query_name)
reader.close()

on peut faire aussi simplement for align in reader:
reader.header : objet représentant le header, de type pysam.libcalignmentfile.AlignmentHeader, sur lequel on peut faire : reader.header.to_dict()
reader.header.get('SQ') : la liste des sequences de la banque (du type {'SN': 'NC_000010.11')
reader.references : la liste de tous les noms des séquences références.

si le fichier est trié et indexé, on peut demander une région :

for align in reader.fetch(contig = 'NC_000018.10', start = 10000000, end = 10100000):
  print(align.to_dict())

Ecriture d'un fichier bam :

writer = pysam.AlignmentFile('outFile.bam', 'wb', header = myHeader) où myHeader est par exemple obtenu d'un autre fichier (myHeader = reader.header)

writer = pysam.AlignmentFile('toto.bam', 'wb', header = reader.header)
  writer.write(align)
writer.close()

objets correspondant à une ligne du bam, de type pysam.libcalignedsegment.AlignedSegment :

align.to_dict() : renvoie un dictionnaire correspondant à la ligne.
align.flag : le flag (2ème colonne du bam)
align.is_mapped : True si la sequence a été mappée.
align.is_forward : True si la séquence a été alignée en forward.
align.is_reverse : True si la séquence a été alignée en reverse.
align.is_read1 et align.is_read2 : True si c'est R1 ou R2 respectivement.
align.is_secondary : True si c'est un alignement secondaire.
align.is_supplementary : True si c'est un alignement supplémentaire.
align.query_name : le nom de la query
align.reference_name : le nom de la référence.
align.query_sequence : la sequence telle qu'elle est dans le bam
align.query_alignment_start, align.query_alignment_start : les positions de la partie alignée, origine à 0, inclus pour le start, exclu pour le end.
align.reference_alignment_start, align.reference_alignment_start : les positions de la partie alignée, origine à 0, inclus pour le start, exclu pour le end.
align.
align.mapping_quality : le MAPQ
align.mate_is_mapped, align.mate_is_forward, align.mate_is_reverse : propriété du mate.
align.get_reference_sequence() : la séquence de la référence alignée comme elle est dans le bam.
align.cigarstring : le CIGAR
align.cigartuples : liste de (operation, position) avec un code pour chaque opération (7 pour =, 1 pour I, 2 pour D, ...
align.get_forward_sequence() : le read original (il peut avoir été reverse-complémenté dans le bam).
align.get_tags() : une liste des (tag, valeur)
align.get_tag('XS') : la valeur associée au tag XS.
align.has_tag('XS') : teste si le tag est présent.

On peut exécuter des commandes samtools :

pysam.sort('-o', 'output.bam', 'input.bam') : tri d'un fichier.
pysam.index('myFile.bam') : indexation d'un fichier trié.

On peut faire un pileup sur une région et itérer sur les positions :

reader = pysam.AlignmentFile('myFile.bam', 'rb')
for pileupcolumn in reader.pileup("chr1", 100, 120):
  ...

programmer en python, tutoriel python, graphes en python, Aymeric Duclert