perl - 如何通过perl中的for循环控件打开和读取文件
问题描述
我在一个文件夹中有很多文件。我想根据参考文件按顺序打开和阅读它们。我的文件名:
AAAAA_AAAAA.CCCCC3.1.bbb.DDDDD.1.fa
AAAAA_AAAAA.CCCCC3.1.bbb.DDDDD.2.fa
AAAAA_AAAAA.CCCCC3.1.bbb.DDDDD.3.fa
AAAAA_AAAAA.CCCCC3.1.bbb.DDDDD.4.fa
.
.
.
参考文件结构:
chr1 744 745
chr1 1208 1209
chr2 1250 1251
chr2 1454 1455
chr3 1676 1677
chr3 1683 1684
输入文件结构:
AAAAA_AAAAA.CCCCC3.1.bbb.DDDDD.1.fa
>1 dna:
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNTATGTGAGAAGATAGCTGAA
CGCCTTGTCCACATCATCTTACTGCTGAGAGTTGAGCTCACCCTCAGTCCCTCACAGTTC
AAAAA_AAAAA.CCCCC3.1.bbb.DDDDD.2.fa
>2 dna:
GAGAGCTGGCTTCTAGGCATGCTTCCTTTTGAGAGCTGAGGACAGGACAGAACCCTCCCG
CATCCTGCCTGACTGTAGACGTACCTGCTAACCTCCTCATGTTAGTGGCTGGGATAGATT
GTGGGAAAAGCATGTGTAAGCATTGGGCCTGAACTCCCGTGTATCTGAGTTGAATACAGC
GATTTCCAACATCCTTCTTCAATAGGAGTGTAGCTAGGTTCCAACTCCCATGTCCGAGTG
GGTAGCAGACATCTGCCTTCCATGCATACACACTTCTGAGAGTTGAGCTTATGGCCTGTA
ACCCTACCTCCTGCCTGCAGCTACCTTTTGCTTCCAAAAGTCCTAGGCTCGCTGCTTCAC
CAAAGTGTTGGGAGAGGTAACTGTTGTCTCCCGGCACACAAGACTAGTGCCTCCAAGCTC
AATCCAGCGATTTCCCAGTAATTCCTGGGTTAGACTGGTGCTACATACTAAGTTCCATAC
GTGAGTAGGTAGTTGAAAGCCTTGTCCAAAAACATCTTACTTCTGAGAGTTGAGCTCACC
CTCAGTCCCTCACAGTTCCACACTGCCTGCAGAGTGAGTTTCCCACGTCTTCATCAGAGA
CTTTTGCCAGAGGCTTCTGAGACGCAAGTTAACAATGCAAACAGGAGGGTATACCCAGGT
GCAGTAGATTGGTTATCTGGGAACCTCCTTACTCAGAATACTGTTACCTTCACACTGTCA
TAAGAATGCAGCTAGTTGAGAGCTGGCTTCTAGGCATGCTTCCCTGTGAGAGCTGAGGAC
我的输出:
chr1 A
chr1 G
chr2 C
chr2 C
chr3 T
chr3 T
我可以使用 bioperl 来查找位置并一一打印出值(逐个文件)。
然后我尝试从文件夹中打开和读取文件。
my $dir = '/home/Documents/Folder/';
opendir(DIR, $dir) or die $!;
my @files = grep (/.fa$/, readdir(DIR));
for my $list(@files){ ##try to get the last number from file name##
my @lines = split /\./, $list}
打开并阅读我的参考文件
open my $POS, '<', 'CanFam3_SNP_POS.txt' or die $!;
我将所有文件放入一个数组并对其进行排序。
my @sorted = @files;
foreach my $i (0..$#sorted)
然后我尝试使用循环控件根据参考文件第 1 列的值打开和读取文件。例如 chr1,应读取并处理 AAAA_AAAAA.CCCCC3.1.bbb.DDDDD.1.fa。如果从参考文件中读取 chr2,则中断循环,然后打开并读取 AAAAA_AAAAA.CCCCC3.1.bbb.DDDDD.2.fa,使用 chr2 处理文件。
open my $fh, '<', "/home/Documents/Folder/$sorted[$i]" or die $!;
while (my $line = <$POS>){
chomp($line);
if ($line =~ /chr$lines[5]/g){
my @positions = split (/\t/, $line);
print "$positions[0]","\t","$positions[1]","\t", substr($so->seq(),
$positions[1], $positions[2] - $positions[1]),"\n";
last if ($line !~ /chr$lines[5]/g)
}
}
我想我对这段代码有一些问题。我可以使用 perl 来完成这个过程吗?我误解了一些观点吗?
解决方案
这里的关键是只查看参考文件中给定行所需的特定 FastA 文件。从您发布的代码摘录来看,您似乎正在尝试读取每一行的每个文件(并且没有设法做到这一点)。
因此,请考虑以下事项:
#!/usr/bin/perl
use warnings;
use strict;
use autodie;
use feature qw/say/;
use File::Basename;
# Map the fasta files in a given directory to chr numbers.
my $fa_dir = '/home/Documents/Folder/';
my %fa_files =
map { (split /\./, fileparse($_, '.fa'))[5] => $_ } glob("$fa_dir/*.fa");
open my $chrs, '<', 'CanFam3_SNP_POS.txt';
# Read each line of the reference file
while (<$chrs>) {
chomp;
# Split up the chr and offsets
my @fields = split /\s+/, $_; #/
# Extract the chr number
my $chr = $fields[0];
$chr =~ s/^chr//; #/
warn "Unknown chr $chr!\n" and next unless exists $fa_files{$chr};
# And read from the appropriate fasta file
# You should probably use a library to read the file instead of
# this mess, but I don't know which ones are good. Based on your code
# you might be trying to use one already?
open my $fa, '<', $fa_files{$chr};
my $hdr = <$fa>;
my $data = join "", <$fa>;
$data =~ s/[^ACGT]+//sg;
close $fa;
# And display the requested part
warn "Invalid offset for chr $chr\n" and next unless length($data) > $fields[1];
my $range = substr $data, $fields[1], $fields[2] - $fields[1];
say "chr$chr $range";
}
它将给定目录中的每个.fa
文件存储到一个哈希表中,由文件名的最后一个元素作为键,该元素对应chr
于参考文件中的内容。这使得查找您需要读取的确切文件以打印出请求的片段变得容易。
还要注意使用glob()来读取文件名,而不是opendir()
/ readdir()
。以这种方式基于扩展名进行过滤更简单,并使用File::Basename以独立于操作系统的方式仅获取文件名减去路径和扩展名。
推荐阅读
- pycairo - 尝试测试 manim 示例时出错
- reactjs - 通过公共 IP 打开时控制台上的 Django 部署错误
- nfc - 使用 DESfire EV1/EV2 进行外部身份验证?
- java - ScrollView 不完全滚动孩子
- javascript - React Storybook 中的 TaskHarness
- ruby-on-rails - Docker 入口点总是返回权限被拒绝
- php - 在 Bootstrap 5 中更改主按钮颜色
- cmake - “通过魔法”将不正确的依赖项添加到 CMake 目标,variable_watch 报告没有写入
- datetime - 离线模式下的otp当前日期时间问题
- api - VS Code API 调用不适用于基本授权