如何从BAM文件中提取fastq

最新推荐文章于 2024-06-25 01:27:02 发布
最新推荐文章于 2024-06-25 01:27:02 发布
阅读量2.4k 收藏 9
点赞数 4
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-NC-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/u012110870/article/details/102804236
版权

虽然高通量测序分析最常用的操作是将fastq比对到参考基因组得到BAM文件，但偶尔我们也需要提取BAM文件中特定区域中fastq。最开始我认为这是一个非常简单的操作，因为samtools其实已经提供了相应的工具samtools fastq.

以biostar handbook的Ebola病毒数据为例，首先获取比对得到的BAM文件。

# 建立文件夹
mkdir -p refs
# 根据Accession下载
ACC=AF086833
REF=refs/$ACC.fa
efetch -db=nuccore -format=fasta -id=$ACC | seqret -filter -sid $ACC > $REF
bwa index $REF
SRR=SRR1553500
fastq-dump -X 100000 --split-files $SRR
BAM=$SRR.bam
R1=${SRR}_1.fastq
R2=${SRR}_2.fastq
TAG="@RG\tID:$SRR\tSM:$SRR\tLB:$SRR"
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM
samtools index $BAM

让我们尝试提取下前1kb比对的fastq文件

samtools view -h SRR1553500.bam KJ660346.1:1-1000 | samtools fastq -1 read_1.fq -2 read_2.fq -
#[M::bam2fq_mainloop] discarded 0 singletons
#[M::bam2fq_mainloop] processed 8702 reads

看起来似乎我们成功了，那让我们把这些序列比对回去吧

bwa mem refs/AF086833.fa read_1.fq read_2.fq
[mem_sam_pe] paired reads have different names: "SRR1553500.11772", "SRR1553500.43775"

你会发现这里产生了错误。这个错误说的是配对的read名字不同。导致问题的原因是：samtools fastq会按照fastq在bam里原本的顺序进行提取，而按照位置排序的bam的两个配对read并不是前后紧挨着的。所以我们需要进行一波排序。

cat read_1.fq \
| paste - - - - \
| sort -k1,1 -S 3G \
| tr '\t' '\n' \
| gzip > read_1.fq.gz
cat read_2.fq \
| paste - - - - \
| sort -k1,1 -S 3G \
| tr '\t' '\n' \
| gzip > read_2.fq.gz

这里将5个shell命令用管道连起来完成这个任务，咋看起来比较复杂，但是理解每一步的目的，每个参数的作用，以后就能活用起来了。

首先用cat逐行读取read，传入到paste中。paste的作用是将多个文件合并成一个文件，做法比较简单粗暴，举个例子

cat A.txt
A
B
C
cat B.txt
D
E
F
paste A.txt B.txt
A   D
B   E
C   F

cat read_1.fq | paste - - - -表示读取4行合并成一行。相当于fastq里的read都会用一行而不是四行表示。接着用sort按照第一行(-k1,1)排序，内存的缓冲数据为3G(-S 3G). 然后将制表符替换成换行符，让一行又变回到4行，后续是压缩成gz，降低磁盘占用。

这样做其实还没又解决问题，因为这两个文件里的reads数其实不相同，我们需要进一步从里面剔除掉只有一方才有的read。这里使用的seqkit，一个非常强大的序列处理轮子。

seqkit grep -f <(seqkit seq -n -i read_1.fq.gz ) read_2.fq.gz  > read_flt_2.fq
seqkit grep -f <(seqkit seq -n -i read_2.fq.gz ) read_1.fq.gz  > read_flt_1.fq
wc -l read_flt_1.fq
   15028 read_flt_1.fq
wc -l read_flt_2.fq
   15028 read_flt_2.fq

终于，你得到了配对存放的两个fastq文件。

通过这个案例我想表达的观点是：虽然有些时候没有现成的工具可以解决你的问题，但是你可以通过组合几个现有的工具来完成。也就是将一个大问题拆解成几个独立的小问题，然后逐个击破。

当然，你首先得熟悉一些已有的工具，比如说常用unix命令:sort, paste, tr等，和一些非常强大生信命令行工具：samtools, bedtools, seqkit等。

一个更佳方便的方法

前面只是问题的一种方法，可能还不是最好的，但是我当初最先想到的一个比较粗暴的方法。当然更佳巧妙的方法是，提取序列之后可以按照read name排序，然后提取

samtools view -h SRR1553500.bam KJ660346.1:1-1000 | samtools sort -n | samtools fastq -1 read_1.fq -2 read_2.fq -s singleton.fq -

是不是更好理解了。

Basic local alignment search tool (BLAST)

包括：blastn, blastp, blastx, tblastn, tblastx等. 使用conda安装即可。

conda install -c bioconda blast
# blast安装perl模块的方法
conda install perl-digest-md5

BLAST的主要理念

  • Search may take place in nucleotide and/or protein space or translated spaces where nucleotides are translated into proteins.
  • Searches may implement search “strategies”: optimizations to a certain task. Different search strategies will return different alignments.
  • Searches use alignments that rely on scoring matrices
  • Searches may be customized with many additional parameters. BLAST has many subtle functions that most users never need.

本地BLAST的基本步骤

  1. 用makeblastdb为BLAST提供数据库
  2. 选择blast工具，blastn,blastp
  3. 运行工具，有必要的还可以对输出结果进行修饰

第一步：建立检索所需数据库

BLAST数据库分为两类，核酸数据库和氨基酸数据库，可以用makeblastbd创建。可以用help参数简单看下说明。

$ makeblastdb -help
USAGE
  makeblastdb [-h] [-help] [-in input_file] [-input_type type]
    -dbtype molecule_type [-title database_title] [-parse_seqids]
    [-hash_index] [-mask_data mask_data_files] [-mask_id mask_algo_ids]
    [-mask_desc mask_algo_descriptions] [-gi_mask]
    [-gi_mask_name gi_based_mask_names] [-out database_name]
    [-max_file_sz number_of_bytes] [-logfile File_Name] [-taxid TaxID]
    [-taxid_map TaxIDMapFile] [-version]
-dbtype <String, `nucl', `prot'>

具体以拟南芥基因组作为案例，介绍使用方法：
注： 拟南芥的基因组可以在TAIR上下在，也可在ensemblgenomes下载。后者还可以下载其他植物的基因组

# 下载基因组
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gzip -d Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# 构建核酸BLAST数据库
makeblastdb -in Arabidopsis_thaliana.TAIR10.dna.toplevel.fa -dbtype nucl -out TAIR10 -parse_seqids

# 下载拟南芥protein数据
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
# 构建蛋白BLAST数据库
gzip -dArabidopsis_thaliana.TAIR10.pep.all.fa.gz
makeblastdb -in  Arabidopsis_thaliana.TAIR10.pep.all.fa -dbtype prot -out TAIR10 -parse_seqids

如果你从NCBI或者其他渠道下载了格式化过的数据库，那么可以用blastdbcmd去检索blast数据库，参数很多，常用就如下几个

  • db string : string表示数据库所在路径
  • dbtype string,: string在（guess, nucl, prot)中选择一个
  • 检索相关参数
    • -entry all 或 555, AC147927 或 gnl|dbname|tag'
    • -entry_batch 提供一个包含多个检索关键字的文件
    • -info 数据库基本信息
  • 输出格式 -outfmt %f %s %a %g ...默认是%f
  • out 输出文件
  • show_blastdb_search_path: blast检索数据库路径

使用案例

# 查看信息
blastdbcmd -db TAIR10 -dbtype nucl -info
# 所有数据
blastdbcmd -db TAIR10 -dbtype nucl -entry all | head
# 具体关键字，如GI号
blastdbcmd -db TAIR10 -dbtype nucl -entry 3 | head

还有其他有意思的参数，可以看帮助文件了解

可选：BLAST安装和更新nr和nt库

安装nt/nr库需要先进行环境变量配置，在家目录下新建一个.ncbirc配置文件，然后添加如下内容

; 开始配置BLAST
[BLAST]
; 声明BLAST数据库安装位置
BLASTDB=/home/xzg/Database/blast
; Specifies the data sources to use for automatic resolution
; for sequence identifiers
DATA_LOADERS=blastdb
; 蛋白序列数据库存本地位置
BLASTDB_PROT_DATA_LOADER=/home/xzg/Database/blast/nr
; 核酸数据库本地存放位置
BLASTDB_NUCL_DATA_LOADER=/home/xzg/Database/blast/nt
[WINDOW_MASKER]
WINDOW_MASKER_PATH=/home/xzg/Database/blast/windowmasker

配置好之后，使用BLAST+自带的update_blastdb.pl脚本下载nr和nt等库文件（不建议下载序列文件，一是因为后者文件更大，二是因为可以从库文件中提取序列blastdbcmd -db nr -dbtype prot -entry all -outfmt "%f" -out nr.fa ，最主要是建库需要花费很长时间），直接运行下列命令即可自动下载。

nohup time update_blastdb.pl nt nr > log &

如果你不像通过update_blastdb.pl下载nr和nt等库文件，也可以是从ncbi上直接下载一系列nt/nr.xx.tar.gz文件，然后解压缩即可，后续还可以用update_blastdb.pl进行数据更新。

报错： 使用update_blastdb.pl更新和下载数据库时候出现模块未安装的问题。解决方法，首先用conda安装对应的模块，然后修改update_blastdb.pl的第一行，即shebang部分，以conda的perl替换，或者按照如下方法执行。

perl `which update_blastdb.pl`

下载过程中请确保网络状态良好，否则会出现Downloading nt.00.tar.gz...Unable to close datastream报错。

第二步：选择blast工具

根据不同的需求，比如说你用的序列是氨基酸还是核苷酸，你要查找的数据是核甘酸还是氨基酸，选择合适的blast工具。不同需求的对应关系可以见下图（来自biostars handbook）

2013053-7eda5a0afe83a170.png
BLAST工具

不同工具的应用范围虽然不同，但是基本参数都是一致的，学会blastn，基本上其他诸如blastp,blastx也都会了。

blastn的使用参数很多 blastn [-h] ，但是比较常用是如下几个

  • -db : 数据库在本地的位置，或者是NCBI上数据库的类型,如 -db nr
2013053-9013c614e8290e5d.png
BLAST库
  • -query: 检索文件
  • -query_loc : 指定检索的位置
  • -strand： 搜索正义链还是反义链，还是都要
  • out : 输出文件
  • -remote： 可以用NCBI的远程数据库， 一般与 -db nr
  • -evalue 科学计数法，比如说1e3，定义期望值阈值。E值表明在随机的情况下，其它序列与目标序列相似度要大于这条显示的序列的可能性。 与S值有关，S值表示两序列的同源性，分值越高表明它们之间相似的程度越大

E值总结： 1.E值适合于有一定长度，而且复杂度不能太低的序列。2. 当E值小于10-5时，表明两序列有较高的同源性，而不是因为计算错误。3. 当E值小于10-6时，表时两序列的同源性非常高，几乎没有必要再做确认。

  • 与联配计分相关参数： -gapopen，打开gap的代价；-gapextend， gap延伸的代价；-penalty：核酸错配的惩罚； -reward， 核酸正确匹配的奖励；

  • 结果过滤：-perc_identity, 根据相似度

注 BLAST还提供一个task参数，感觉很有用的样子，好像会针对任务进行优化速度。

第三步：运行blast，调整输出格式。

我随机找了一段序列进行检索

echo '>test' > query.fa
echo TGAAAGCAAGAAGAGCGTTTGGTGGTTTCTTAACAAATCATTGCAACTCCACAAGGCGCCTGTAATAGACAGCTTGTGCATGGAACTTGGTCCACAGTGCCCTACCACTGATGATGTTGATATCGGAAAGTGGGTTGCAAAAGCTGTTGATTGTTTGGTGATGACGCTAACAATCAAGCTCCTCTGGT >> query.fa

用的是blastn 检索核酸数据库。最简单的用法就是提供数据库所在位置和需要检索的序列文件

blastn -db BLAST/TAIR10 -query query.fa
# 还可以指定检索序列的位置
blastn -db BLAST/TAIR10 -query query.fa  -query_loc 20-100
# 或者使用远程数据库
blastn -db nr -remote -query query.fa
blastn -db nt -remote -query query.fa

以上是默认输出，blast的-outfmt选项提供个性化的选择。一共有18个选择，默认是0。
0 = Pairwise,
1 = Query-anchored showing identities,
2 = Query-anchored no identities,
3 = Flat query-anchored showing identities,
4 = Flat query-anchored no identities,
5 = BLAST XML,
6 = Tabular,
7 = Tabular with comment lines,
8 = Seqalign (Text ASN.1),
9 = Seqalign (Binary ASN.1),
10 = Comma-separated values,
11 = BLAST archive (ASN.1),
12 = Seqalign (JSON),
13 = Multiple-file BLAST JSON,
14 = Multiple-file BLAST XML2,
15 = Single-file BLAST JSON,
16 = Single-file BLAST XML2,
17 = Sequence Alignment/Map (SAM),
18 = Organism Report
其中6，7，10，17可以自定输出格式。默认是

qaccver saccver pident length mismatch gapopen qstart qend sstart send evalue bitscore

简写含义
qaccver查询的AC版本（与此类似的还有qseqid,qgi,qacc,与序列命名有关）
saccver目标的AC版本（于此类似的还有sseqid,sallseqid,sgi,sacc,sallacc,也是序列命名相关)
pident完全匹配百分比 （响应的nident则是匹配数）
length联配长度(另外slen表示查询序列总长度，qlen表示目标序列总长度)
mismatch错配数目
gapopengap的数目
qstart查询序列起始
qstart查询序列结束
sstart目标序列起始
send目标序列结束
evalue期望值
bitscoreBit得分
score原始得分
AC：accession

以格式7为实例进行输出，并且对在线数据库进行检索

blastn -task blastn -remote -db nr -query query.fa  -outfmt 7 -out query.txt
head -n 15 query.txt
2013053-beedfae7aaa0cd45.png
BLAST BEST HIT

如果想输入序列，增加对应的格式qseq, sseq

blastn -query query.fa -db nr -outfmt ' 6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send qseq sseq evalue bitscore'

对已有序列进行注释时常见的best hit only模式命令行

blastn -query gene.fa -out gene.blast.txt -task megablast -db nt -num_threads 12 -evalue 1e-10 -best_hit_score_edge 0.05 -best_hit_overhang 0.25 -outfmt "7 std stitle" -perc_identity 50 -max_target_seqs 1
# 参数详解
-task megablast : 任务执行模式，可选有'blastn' 'blastn-short' 'dc-megablast' 'megablast' 'rmblastn'
-best_hit_score_edge 0.05 :  Best Hit 算法的边界值，取值范围为0到0.5，系统推荐0.1
-best_hit_overhang 0.25 : Best Hit 算法的阈值，取值范围为0到0.5，系统推荐0.1
-perc_identity 50 : 相似度大于50
-max_target_seqs 1 : 最多保留多少个联配

仅仅看参数依旧无用，还需要知道BLAST的Best-Hits的过滤算法。假设一个序列存在两个match结果,A和B，无论A还是B，他们的HSP（High-scoring Segment Pair, 没有gap时的最高联配得分）一定要高于best_hit_overhang，否则被过滤。如果满足下列条件则保留A

  • evalue(A) >= evalue (B)
  • score(A)/length(A) < (1.0-score_edge)*score(B)/length(B)

搭建网页BLAST

曾经的BLAST安装后提供wwwblast用于构建本地的BLAST网页工具，但是BLAST+没有提供这个工具，好在BLAST足够出名，也就有人给它开发网页版工具。如viroBLAST和Sequenceserver, 目前来看似乎后者更受人欢迎。

有root安装起来非常的容易

sudo apt install ruby ncbi-blast+ ruby-dev rubygems-integration npm
sudo gem install sequenceserver

数据库准备,前面步骤已经下载了拟南芥基因组的FASTA格式数据

sequenceserver -d /到/之前/建立/BLAST/文件路径

然后就可以打开浏览器输入IP:端口号使用了。

Sequenceserver高级用法

开机启动：

新建一个/etc/systemd/system/sequenceserver.service文件，添加如下内容。注意修改ExecStart.

[Unit]
Description=SequenceServer server daemon
Documentation="file://sequenceserver --help" "http://sequenceserver.com/doc"
After=network.target

[Service]
Type=simple
User=seqservuser
ExecStart=/path/to/bin/sequenceserver -c /path/to/sequenceserver.conf
KillMode=process
Restart=on-failure
RestartSec=42s
RestartPreventExitStatus=255

[Install]
WantedBy=multi-user.target

然后重新加载systemctl

## let systemd know about changed files
sudo systemctl daemon-reload
## enable service for automatic start on boot
systemctl enable sequenceserver.service
## start service immediately
systemctl start sequenceserver.service

nginx反向代理：我承认没有基本的nginx的知识根本搞不定这一步，所以我建议组内使用就不要折腾这个。简单的说就是在nginx的配置文件的server部分添加如下内容即可。

location /sequenceserver/ {
    root /home/priyam/sequenceserver/public/dist;
    proxy_pass http://localhost:4567/;
    proxy_intercept_errors on;
    proxy_connect_timeout 8;
    proxy_read_timeout 180;
}

假如你有docker的话

# ubuntu
docker run --rm -itp 4567:4567 -v /path-to-database-dir:/db wurmlab/sequenceserver:1.0.11
# centos
docker run --privileged --rm -itp 4567:4567 -v /path-to-database-dir:/db wurmlab/sequenceserver:1.0.11

注意，你得在宿主机器上/path-to-database-dir建好blast索引，或者用docker exec -it docker容器名 /bin/bash 进入到容器里，用sequenceserver -d db &运行

参考资料：http://www.sequenceserver.com/doc/

术语列表

引自BLAST Glossary

参考资料

本文还可在http://xuzhougeng.top/阅读

版权声明：本博客所有文章除特别声明外，均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。

2013053-de4f1abd858af2a1.png
扫码即刻交流
徐洲更hoptop
关注
生信分析进阶4 - 比对结果的FLAG和CIGAR信息含义与BAM文件指定区域提取
LittleComputerRobot的博客
06-01 377
BAM比对文件FLAG和CIGAR含义，samtools提取指定BAM区域。
生信技能57 - Samtools获取指定外显子区域depth和提取BAM文件序列
最新发布
LittleComputerRobot的博客
08-24 106
Samtools获取指定外显子区域depth和提取BAM文件序列
bam2fastq:从bam到FASTQ的简单转换器
05-24
bam2fastq bam2fastq是一个从BAM文件中提取序列和质量的程序。 原始版本可以在找到。 随后对其进行了修改，以处理具有成对和不成对读和写到stdout混合的BAM文件。 安装 需要make，gcc和zlib压缩库-所有这些都应该存在于大多数类Unix系统（包括Mac）上。 首先从github克隆存储库，包括其依赖项： git clone --recursive https://github.com/jts/bam2fastq 然后转到目录并运行make： cd bam2fastq make 作者 该程序最初由Phillip Dexheimer在HudsonAlpha开发。 它打包在github上，并由Jared Simpson稍加修改。
如何高效地从BAM文件中提取fastq
xuzhougeng blog
04-21 1956
在一年前，我写过一篇文章，叫做如何从BAM文件中提取fastq，之前也发现了从BAM里面提取Fastq是有些麻烦，只不过最后通过samtools的子命令实现了数据提取，实现功能之后也没有再去思考如何提高效率。 最近读到每周文献-190419-植物单细胞BAM重比对以及假基因研究时，发现里面提到了一个工具叫做 bazam, 功能就是提取Fastq文件，文章发表在 Genome Bio...
bwa shm 测试
u014210048的专栏
06-25 133
【代码】bwa shm 测试。
格式转换：BAM 转 FASTQ
鬼话
04-03 5579
有时候，分析已发表数据的时候，避免不了会遇到作者上传的数据是BAM格式，但是作者用的基因组又不是我想要的，所以我就需要将BAM转换为FASTQ，加上许多分析工具都是需要以FASTQ文件为起始输入，据说BAM转换为FASTQ是生信中的一个常见步骤，虽然现在还没遇到，提前整理，以备不时之需😑工具有很多，接下来就罗列几个工具。
bam格式转换为Fastq/Fasta格式
热门推荐
lazymark2的博客
09-18 1万+
bam格式转换为Fastq/Fasta格式 Samtools Fastq GATK SamToFastq Bedtools bamtofastq 举例说明，比如说我们现在有一个转录组比对文件D1_D1.sort.bam： samtools view -H D1_D1.sort.bam | tail @SQ SN:19 LN:58617616 @SQ SN:20 LN:64444167 @SQ SN:21 LN:46709983 @SQ SN:22 LN:50818468 @SQ SN:X LN:1560
bam获取序列_如何从BAM文件中提取fastq
weixin_39722196的博客
12-20 2202
虽然高通量测序分析最常用的操作是将fastq比对到参考基因组得到BAM文件，但偶尔我们也需要提取BAM文件中特定区域中fastq。最开始我认为这是一个非常简单的操作，因为samtools其实已经提供了相应的工具samtools fastq.以biostar handbook的Ebola病毒数据为例，首先获取比对得到的BAM文件。# 建立文件夹mkdir -p refs# 根据Accession下载...
linux怎么根据比对结果的bam文件按照所要的位置区间进行reads的提取成fastq的格式，具体代码流程是什么，或者是已经提成bam再怎么把截取的bam转成fastq
06-09
1. 使用samtools view命令根据指定的位置区间将BAM文件中的reads提取出来，并将结果保存为一个新的BAM文件。 ``` samtools view -b input.bam chr:start-end > output.bam ``` 其中，chr表示染色体名称，start和...
bam文件读取_bam格式文件处理大全(一）
weixin_29422697的博客
12-30 3886
sam文件是短序列比对生成的文件，是二代测序中最核心的文件。在RNAseq，变异检测等分析中，都需要首先生成sam文件格式。bam文件是sam格式的二进制格式，转换为二进制之后，可以减小文件的存储。掌握sam/bam文件的操作是处理二代测序数据的非常重要的内容，例如sam与bam的转换，排序，建立索引，reads计数等等操作。1 生成sam文件有很多软件可以生成sam格式，例如bwa，bowie2...
bam文件转fq.gz文件
S_AGZX的博客
06-19 2870
fastq文件格式；bam文件格式；bam转换fastq的方法；代码示例
踩坑日记|bam文件转为fastq文件
XIAOWANGaxs的博客
12-10 1135
好了，今天的踩坑日记就到这里了，其它实现 bam2fastq 工具（例如：Picard、GATK、Biopython）没有做过尝试，大家尝试一下发现的坑，可以留言交流一下。命令进行 bam 文件转换为 fastq 文件，可以对单端比对的生成的 bam 和双端比对生成的 bam 进行转换，详细使用方法和参数见。命令进行 bam 文件转换为 fastq 文件，并且可以将单端比对、双端比对的 bam 文件转换为 fastq 文件，使用方法和参数见。使用的一些工具和工具使用过程中发现的问题。
10x单细胞转录组数据bam格式转为fastq格式
weixin_47890536的博客
04-19 2040
在下载NCBI等公共数据库中的单细胞测序数据，发现仅提供了处理后的bam文件时，可以用来将bam转为fastq文件进行重新分析，实现不同数据集上游分析的统一，比如参考基因组的一致等。
unmapbam to fastq和自己的annovar格式~~~
浮生终有醒
05-31 1676
#!perl use warnings; use strict; die "perl $0 \n" if @ARGV != 2; my %hash; open BAM, "samtools view $ARGV[0] |" or die $!; while() { chomp; my @tmp = split; push @{$hash{$tmp[0]}}, "$tmp[1]\t$tm
批量从bam文件获取指定位置的碱基
Cassiel60的博客
10-11 517
批量从bam文件获取指定位置(RS.xlsx)的碱基后输出到一个excel，按照编号进行排序。
生信人迷惑的一天 bam转fq
qq_43337286的博客
07-22 1568
bam格式转fastq的二三事
method by id是什么_根据id快速提取fastq序列
weixin_32464651的博客
12-20 268
喜欢就点关注吧!根据fastq序列的id，从原始fastq中提取序列这个操作，应该是大家在处理序列文件的过程中经常遇到的。如果大家用过Biopython，应该知道Bio模块在做fastq这些文件的处理时非常方便。但是有时序列达到几百万几千万条的时候，Bio的速度可能就无法满足要求了。还是举个例子比较好，我从比对筛选过滤之后的bam文件中提取了第一列序列名，保存为id.name文件，想根据...
Adam学习13之Fasta/Fastq/SAM/BAM文件格式数据读取
Keep Learning
04-30 5702
0.代码（读取方法）： package org.bdgenomics.adamLocal.algorithms.test import org.apache.spark.SparkConf import org.apache.spark.SparkContext import org.apache.spark.sql.SQLContext import org.bdgenomics.adam.
从比对的BAM文件中提取所有突变
04-04
从比对的BAM文件中提取所有突变需要进行以下步骤： 1. 使用samtools将BAM文件转换为SAM文件： ``` samtools view -h input.bam > output.sam ``` 2. 使用samtools将SAM文件转换为BAM文件，并按照坐标排序： ``` ...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值