Seqtk(fq转换为fa，格式化序列，截取序列，随机抽取序列）

下载安装

apt-get install seqtk

seq 序列常规转换

1.将fastq转换成fasta：

seqtk seq -a Sample_R1.fq.gz > Sample_R1.fa

2.将fastq序列做反向互补分析：

seqtk seq -r Sample_R1.fq.gz > Sample_Revc_R1.fq

3.sample 随机抽样

seqtk sample -s100 Sample_R1.fq.gz 10000
4.subseq 提取序列
4.1根据输入的bed文件信息，将固定区域的序列提取出来：

seqtk subseq in.fa reg.bed > out.fa

4.2根据输入的name list，提取相应名称序列：

seqtk subseq in.fq name.lst > out.fq

5.截取序列

5.1切除reads的前5bp，以及后10bp：

seqtk trimfq -b 5 -e 10 in.fq > out.fq

源自生物信息常用软件 seqtk - 知乎 (zhihu.com)

例1，从一个fastq文件里提取10000reads

文件目录/位置 |seqtk sample -s 60 - 500 >100000reads.fq（输出文件）
用

less 100000reads.fq |wc -l
确认一下

显示2000

上面提取的时候是提取500reads为一行，进入文件之后有2000行，确认无误。
less 100000reads.fq 查看文件具体内容。
更多》》》seqtk的使用说明 - 简书 (jianshu.com)

Trimmomatic（测序数据质控）

Trimmomatic 用于去除 Illumina平台的FASTQ序列中的Adapter，根据碱基质量值修整FASTQ序列文件
下载安装

wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.38.zip

解压

unzip Trimmomatic-0.38.zip
测试运行
java -jar ~/Biosofts/Trimmomatic-0.38/trimmomatic-0.38.jar

结果显示

简单测试一下

ILLUMINACLIP：从reads中剪切adapter和其他Illumina特定序列。
SLIDINGWINDOW：执行滑动窗口修剪，一旦窗口内的平均质量低于阈值，则切割。
LEADING：如果低于阈值质量，则在reads起始处剪切碱基
TRAILING：如果低于阈值质量，则在reads末尾处剪切碱基
CROP：将reads从末尾切割为指定长度
HEADCROP：从reads剪切后低于指定长度，则删除
MINLEN：如果reads低于指定长度，则删除
TOPHRED33：将质量得分转换为Phred-33
TOPHRED64：将质量得分转换为Phred-64

java -jar ~/Biosofts/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 /disk1/shares/Seqs/test_7942raw_1.fq.gz  /disk1/shares/Seqs//test_7942raw_2.fq.gz ./trim_out/output_forward_paired.fq.gz ./trim_out/output_forward_unpaired.fq.gz ./trim_out/output_reverse_paired.fq.gz ./trim_out/output_reverse_unpaired.fq.gz ILLUMINACLIP:/disk1/202031107010230/Biosofts/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 MINLEN:75

结果显示

进入disk1/202031107010230/Biosofts

示例显示四个输出文件

参考生信软件 | Trimmomatic (测序数据质控) - 知乎 (zhihu.com)

最后我们使用fastqc查看一下

fastqc结果

Multiqc整合多个质控结果

四个输出文件一起对比

好了，直接下班！！！

2022.928

补充一个列子，用Seqtk提取30个蛋白质氨基酸序列

提取DNA

-寻找一个AA.gff文件,进行如下操作

grep $'\tProtein' GCA_000817325.1_ASM81732v1_genomic.gff|grep -v "pseudogene" |awk -v FS="\t" -v OFS="\t" '{print $1,$4,$5,$7,$9}'|sed 's/\tID.*;locus_tag=/\t/g'|sed 's/;.*;protein_id=/\t/g'|sed 's/;.*$//g'|awk -v FS='\t' -v OFS='\t' '{print $1,$2-1,$3,$5,"0",$4,$6}'>genome.bed

运行效果图1

seqtk subseq /文件目录/AA.fna genome.bed >cds_A1.fna
seqtk cds_A1.fna|less

提取目的数量的蛋白质氨基酸序列

cut -f 7 genome.bed |head -n30>pro_name.list #提取出蛋白质前体名称用与.faa文件对比抽取
seqtk subseq /文件目录/AA.faa pro_name.list >A3.faa
less -S A3.faa

图2