我正在编写一个脚本来分析一些数据。我有几个文件子集,我想计算这些文件中的行数并将结果写入 csv 文件。我将尝试做一个例子。我有这两个文件子集:
sample1.ext
sample1.ext2
sample1.ext3
sample2.ext
sample2.ext2
sample2.ext3
我想计算、和中所有文件包含的行数*.ext
,*.ext2
并将*.ext3
结果写入如下所示的 csv 文件中:
count(sample1.ext), count(sample1.ext2), count(sample1.ext3)
count(sample2.ext), count(sample2.ext2), count(sample2.ext3)
在对 中的第一个文件系列进行计数后*.ext
,我将结果输出到 csv 文件的第一列。如何将 中的第二个计数系列的输出写入*.ext2
同一个 csv 文件的第二列?第三列也一样吗?
感谢大家的回答,我试图将它们改编到我的文件中,但不幸的是我做不到。我发布的示例只是一个例子,我使用数字代替奇怪的扩展名,以便最容易理解问题。你们都明白了,但你们太关注现实中不存在的数字了。我将使用真实文件再次向你们解释。这些文件来自基因组数据到参考基因组的映射。我处理这些数据以清理它们,所以我有三个步骤,其中行数会发生变化。所以文件是:
name.sort.bam
name.mapped.bam
name.rmdup.bam
othername.sort.bam
othername.mapped.bam
othername.rmdup.bam
扩展名 bam 是一个压缩文件。要计算此文件中的行数,有一个特殊的命令行:
samtools view -c (file)
我发现的唯一方法是迭代每个*sort.bam
,*mapped.bam
并*rmdup.bam
为每个写一个 txt 输出,然后将它们粘贴到 csv 文件的末尾。有没有办法避免这三个循环并一起做所有事情?抱歉造成误解,你们都有很好的想法!
答案1
您可以使用这个 Perl 脚本:
#! /usr/bin/perl
use strict;
use warnings;
my @names;
my @files;
@ARGV == 1 || die();
opendir(my $dir, $ARGV[0]) || die $!;
while(readdir($dir)) {
if($_ =~ /(.*)\.(sort|mapped|rmdup)\.bam$/) {
grep(/^$1$/, @names) == 0 && push(@names, $1);
}
}
close($dir);
foreach my $name (sort(@names)) {
my @fields;
push(@fields, $name);
foreach my $extension ("sort", "mapped", "rmdup") {
if(! -f "$ARGV[0]/$name.$extension.bam") {
push(@fields, 0);
print STDERR "'$ARGV[0]/$name.$extension.bam' missing\n";
next;
}
my $count = `<"$ARGV[0]/$name.$extension.bam" wc -l`;
chomp($count);
push(@fields, $count)
}
print(join(", ", @fields)."\n")
}
将其保存在系统中的某个位置,使其可执行,并将目录作为参数传递来运行它:
path/to/script path/to/directory
% tree directory
directory
├── name.mapped.bam
├── name.rmdup.bam
├── name.sort.bam
├── othername.mapped.bam
├── othername.rmdup.bam
└── othername.sort.bam
0 directories, 6 files
% perl script.pl directory
name, 0, 0, 0
othername, 0, 0, 0
% for f in directory/*.sort.bam; do printf 'line\n' >>"$f"; done
% perl script.pl directory
name, 1, 0, 0
othername, 1, 0, 0
该脚本的作用是:
- 遍历中的所有文件
path/to/directory
;如果文件名与匹配.*\.(sort|mapped|rmdup)\.bam$
,则将字符串附加在前面.sort.bam
,如果尚未在列表中,.mapped.bam
则将.rmdup.bam
其附加到列表中;@names
- 对于排序
@names
列表中的每个名称$name
,将其附加$name
到列表@fields
;对于中的每个扩展名sort
,mapped
以及rmdup
检查$extension
是否$name.$extension.bam
存在于中path/to/directory
;如果文件不存在,则附加0
到@fields
,打印错误消息并转到下一个$extension
/$name
;如果文件存在,则将 的输出附加<"$name.$extension.bam" wc -l
到@fields
;一旦所有可能的值都已迭代,则打印一行包含已加入$extension
的元素。@fields
,
答案2
假设您想要42, 19, 10207, 3
在每一行上输出类似内容(没有文件名),wc
并且一些Bash
ing 将解决您的问题。
outfile="Result.csv"
for samplenum in $( seq 1 100 ) ; do
line=""
for file in sample${samplenum}.* ; do
numlines=$( wc -l <$file )
line="$line $numlines,"
done
# remove the final comma
line=${line%,}
# not quoting $line below will suppress the initial blank
echo $line >> $outfile
done
阅读man bash
,,然后重读man wc
man seq
man bash
回复评论:
您读过这些man
页面吗?
$( seq 1 100)
被命令的结果所取代seq 1 100
,该命令仅输出从 1 到 100 的整数(阅读man seq
会告诉您)。将其替换为提供您拥有的样本数量的内容。
将代码放入文件(例如test.sh
)中,然后使用运行它以bash -x test.sh
查看详细信息。将测试替换seq 1 100
为seq 1 2
,以避免输出大量内容。
samplenum
保存样本编号,在此示例中,从 1 到 100。
sample
中sample${samplenum}.*
只是一个字符串。它与的值samplenum
和字符串连接起来.*
,产生文件名模式,例如sample1.*
第一次循环for samplenum ...
,sample2.*
第二次循环,等等。
您是否已阅读并理解、、、man bash
然后man wc
再读man seq
一遍man bash
?
答案3
python 选项
有趣的问题。这是应用 Python 的好时机groupby()
由于您的文件位于单个“平面”目录中:
#!/usr/bin/env python3
from itertools import groupby
import os
import sys
dr = sys.argv[1]
# list the files in the directory, split into "sortable" elements
flist = [[item, item.split(".", 1)] for item in os.listdir(dr)]
# sort the file list by first section (until the first found dot)
flist.sort(key=lambda x: x[1][0])
# create sub groups of the files, grouped by first section of name
for key, line in groupby(flist, lambda x: x[1][0]):
line = list(line)
# sort the files by second section of name for correct order in the csv lines
line.sort(key=lambda x: x[1][1])
# count the lines of the files, arrange the csv file
print((", ").join([str(len(open(dr+"/"+f[0]).readlines())) for f in line]))
怎么运行的
如果目录包含九个文件:
sample1.ext 2 lines
sample1.ext2 3 lines
sample1.ext3 3 lines
sample2.ext 1 lines
sample2.ext2 1 lines
sample2.ext3 4 lines
sample3.ext 6 lines
sample3.ext2 1 lines
sample3.ext3 4 lines
该脚本列出文件,将每个名称分成两部分,例如:
sample2
和
ext2
由于行的顺序和文件的长度里面这些线条取决于这两个部分的准确排序。
- 然后脚本排序按名称的第一部分对文件进行排序,因为每个文件(具有相似的第一个名称部分)的长度应该是按第一部分分组合并为一行;,,
sample1
等等。sample2
sample3
csv
随后,创建子组(每行),并按以下方式正确排序:第二名称部分,使(行)号在行中以正确的顺序出现
结果:
python3 '/home/jacob/Bureaublad/create_csv.py' '/home/jacob/Bureaublad/samples'
2, 3, 3
1, 1, 4
6, 1, 4
如何使用
- 将脚本复制到一个空文件中,另存为
create_csv.py
使用包含你的文件的目录作为参数来运行它
python3 /path/to/create_csv.py /path/to/directory_with_files
重要的提示
只要文件不是巨大的. 如果文件是,采用另一种方法来计算行数将会产生更好的性能。
编辑
根据最新信息,在您的问题中添加了脚本的编辑版本。巧合的是,不需要做太多更改:脚本已经使用以下命令,按找到的第一个点拆分文件名:
item.split(".", 1)
由于名称的最后一部分是.bam
,它等于全部文件,对于排序来说毫无意义。
然后我们只需要替换“旧”的方法来计算文件的行数:
str(len(open(dr+"/"+f[0]).readlines()))
通过(python 实现并集成到脚本中)您提供的命令:
str(subprocess.check_output(["samtools", "view", "-c", dr+"/"+f[0]]).decode("utf-8").strip())
编辑后的脚本
#!/usr/bin/env python3
from itertools import groupby
import os
import sys
import subprocess
dr = sys.argv[1]
# list the files in the directory, split into "sortable" elements
flist = [[item, item.split(".", 1)] for item in os.listdir(dr)]
# sort the file list by first section (until the first found dot)
flist.sort(key=lambda x: x[1][0])
# create sub groups of the files, grouped by first section of name
for key, line in groupby(flist, lambda x: x[1][0]):
line = list(line)
# sort the files by second section of name for correct order in the csv lines
line.sort(key=lambda x: x[1][1])
# count the lines of the files, arrange the csv file
print((", ").join([
str(subprocess.check_output(["samtools", "view", "-c", dr+"/"+f[0]]).decode("utf-8").strip())
for f in line]))
笔记
请注意,行内数字的顺序由名称第二部分的排序顺序决定,例如
mapped.bam, rmdup.bam, sort.bam