用于识别样本中最小染色体区域的Python

Question

我有多个示例文件（> 20），如下所示：

chr startpos    endpos
1   14930   818094
1   818161  31595422
2   35593931    35865807
2   35868158    104785784

我想输出样本中常见的区域。例如。如果样本1有：

1   14900   818000

样本2：

1   15000   605000

样本3：

1   25000   705000

我想输出：

1   25000   605000

我还想包括一个多数规则，例如，如果总共20个样本中的10个具有最小区域 - >输出该区域。即我想让它变得灵活多少需要将区域打印到输出的样本。

有没有人有这个python解决方案？

Answer 1

不确定这是否是针对Unix和Linux stackexchange的问题。这听起来更像是一般编程问题。

但是，我鼓励你研究使用pandas。

您可以将样本文件作为数据框导入，指定制表符描述如下：

import pandas as pd
df = pd.read_csv('/tmp/samplefile.csv',sep='\t')

如果您知道startpos将始终小于endpos，您可以通过获取df['startpos']的最大值和df['endpos']的最小值来找到您正在寻找的输出。

Answer 2

有专门用于处理BED格式数据的工具，您可以通过学习获得良好的服务。 bedtools可能是最常见和最容易拾取的，如果你绝对必须使用python它有一个python wrapper。

multiIntersectBed工具是您可能想要的，并且应该非常简单易用：

用法示例：

== Input files: ==

 $ cat a.bed
 chr1  6   12
 chr1  10  20
 chr1  22  27
 chr1  24  30

 $ cat b.bed
 chr1  12  32
 chr1  14  30

 $ cat c.bed
 chr1  8   15
 chr1  10  14
 chr1  32  34

 $ cat sizes.txt
 chr1  5000

== Multi-intersect the files: ==

 $ multiIntersectBed -i a.bed b.bed c.bed
chr1    6       8       1       1       1       0       0
chr1    8       12      2       1,3     1       0       1
chr1    12      15      3       1,2,3   1       1       1
chr1    15      20      2       1,2     1       1       0
chr1    20      22      1       2       0       1       0
chr1    22      30      2       1,2     1       1       0
chr1    30      32      1       2       0       1       0
chr1    32      34      1       3       0       0       1

== Multi-intersect the files, with a header line (titles are the file names): ==

 $ multiIntersectBed -header -i a.bed b.bed c.bed
 chrom  start   end     num     list    a.bed   b.bed   c.bed
 chr1   6       8       1       1       1       0       0
 chr1   8       12      2       1,3     1       0       1
 chr1   12      15      3       1,2,3   1       1       1
 chr1   15      20      2       1,2     1       1       0
 chr1   20      22      1       2       0       1       0
 chr1   22      30      2       1,2     1       1       0
 chr1   30      32      1       2       0       1       0
 chr1   32      34      1       3       0       0       1

== Multi-intersect the files, with a header line and custom names: ==

 $ multiIntersectBed -header -i a.bed b.bed c.bed -names A B C
 chrom  start   end     num     list    A       B       C
 chr1   6       8       1       A       1       0       0
 chr1   8       12      2       A,C     1       0       1
 chr1   12      15      3       A,B,C   1       1       1
 chr1   15      20      2       A,B     1       1       0
 chr1   20      22      1       B       0       1       0
 chr1   22      30      2       A,B     1       1       0
 chr1   30      32      1       B       0       1       0
 chr1   32      34      1       C       0       0       1

== Multi-intersect the files, showing empty regions (note, requires -g): ==

 $ multiIntersectBed -header -i a.bed b.bed c.bed -names A B C -empty -g sizes.txt
 chrom  start   end     num     list    A       B       C
 chr1   0       6       0       none    0       0       0
 chr1   6       8       1       A       1       0       0
 chr1   8       12      2       A,C     1       0       1
 chr1   12      15      3       A,B,C   1       1       1
 chr1   15      20      2       A,B     1       1       0
 chr1   20      22      1       B       0       1       0
 chr1   22      30      2       A,B     1       1       0
 chr1   30      32      1       B       0       1       0
 chr1   32      34      1       C       0       0       1
 chr1   34      5000    0       none    0       0       0

然后，您可以按第4列过滤，以满足您所规定的区域条件。