使用两个软件对fq数据进行比对，得到比对文件sam/bam，并探索比对结果。
- 三个常用参考基因组数据库：
Ensembl：www.ensembl.org
NCBI：
https://www.ncbi.nlm.nih.gov/projects/genome/gu
ide/human/index.shtml
UCSC：http://www.genome.ucsc.edu/
-
Ensembl官网
image.png
image.png
image.png
这种ftp类型的地址是可以直接用weget命令进行下载的
image.png
## 参考基因组准备:注意参考基因组版本信息
# 下载，Ensembl：http://asia.ensembl.org/index.html
# ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/dna/
# 进入到参考基因组目录
cd /teach/database/genome/Ensembl/Homo_sapiens/GRCh38_release95
# 下载基因组序列
wget ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
——————————————————数据太大就不下载了——————————使用现成的————————————————————————————-
软链接下载好的参考基因组到自己的目录下
-
参考基因组注释文件
每一行有九列
image.png
image.png
image.png
# 下载基因组注释文件
wget ftp://ftp.ensembl.org/pub/release-95/gtf/homo_sapiens/Homo_sapiens.GRCh38.95.gtf.gz
习题
1.fastq与fasta文件转换
应用：加深对两种文件格式的理解，并且fa有后续应用。
#进入fastq目录
less -S SRR1039510_1.fastq.gz | paste - - - - |awk '{print $1"\n"$4}' >SRR1039510_1.fa
#将@替换为>
sed -i 's/@/>/g' SRR1039510_1.fa
image.png
2.使用fa文件做NT比对，得到测序样本的NT比对信息
应用：随机收取2000reads做NT比对，可以发现测序样本是否含有其他物种污染等。
3.从gff或者gft文件中获取基因的ID与symbol对应关系，以及biotype类型
应用：ID与symbol转换本地化，不依赖于第三方工具和软件包，并可以根据biotype类型区分mRNA，lncRNA以及miRNA等信息。
hisat2比对
image.png
image.png
# 进入参考基因组目录
cd ~/database/genome/Ensembl/Homo_sapiens/GRCh38_release95
# Hisat2构建索引
hisat2-build Homo_sapiens.GRCh38_release95.genome.fa Homo_sapiens.GRCh38_release95.genome
# 输入输出定义文件夹
index=~/database/genome/Ensembl/Homo_sapiens/GRCh38_release95/Homo_sapiens.GRCh38_release95.genome.
inputdir=~/project/Human_16-Asthma-Trans/data/cleandata/trim_galore
outdir=~/project/Human_16-Asthma-Trans/Mapping/Hisat2
image.png
# 单个样本比对
hisat2 -p 3 -x ${index} -1 ${inputdir}/SRR1039510_1_val_1.fq.gz -2 ${inputdir}/SRR1039510_2_val_2.fq.gz -S ${outdir}/SRR1039510.Hisat_aln.sam
image.png
第一行：fastq有多少reads，双端的reads：正端和反端各测了一遍
image.png
image.png
# sam转bam
samtools sort -@ 3 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam
# 对bam建索引
samtools index SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sorted.bam.bai
# 多个样本批量进行比对，排序，建索引
# Hisat.sh内容
cat /teach/project/Human-16-Asthma-Trans/data/rawdata/sra/sampleId.txt | while read id
do
echo "hisat2 -p 3 -x ${index} -1 ${inputdir}/${id}_1_val_1.fq.gz -2 ${inputdir}/${id}_2_val_2.fq.gz 2>${id}.log | samtools sort -@ 5 -o ${outdir}/${id}.Hisat_aln.sorted.bam - && samtools index ${outdir}/${id}.Hisat_aln.sorted.bam ${outdir}/${id}.Hisat_aln.sorted.bam.bai"
done >Hisat.sh
# 提交后台运行
nohup sh Hisat.sh >Hisat.log &
# 统计比对情况
multiqc -o ./ SRR*log
subjunc比对
image.png
# 进入参考基因组目录
cd /teach/database/genome/Ensembl/Homo_sapiens/GRCh38_release95
# subjunc构建索引
subread-buildindex -o Homo_sapiens.GRCh38_release95.genome Homo_sapiens.GRCh38_release95.genome.fa
# 输入输出定义文件夹
index=/teach/database/genome/Ensembl/Homo_sapiens/GRCh38_release95/Homo_sapiens.GRCh38_release95.genome
inputdir=/trainee2/Oct12/project/project_test/data/cleandata/trim_galore
outdir=/trainee2/Oct12/project/project_test/Mapping/subjunc
# subjunc单样本比对
subjunc -T 3 -i ${index} -r ${inputdir}/SRR1039510_1_val_1.fq.gz -R ${inputdir}/SRR1039510_2_val_2.fq.gz -o ${outdir}/SRR1039510.Subjunc.bam 1>${outdir}/SRR1039510.Subjunc.log 2>&1
# subjunc多样本比对
cat /teach/data/airway/sra/sampleId.txt | while read id
do
echo "subjunc -T 1 -i ${index} -r ${inputdir}/${id}_1_val_1.fq.gz -R ${inputdir}/${id}_2_val_2.fq.gz -o ${outdir}/${id}.Subjunc.bam 1>${outdir}/${id}.Subjunc.log 2>&1 && samtools sort -@ 6 -o ${outdir}/${id}.Subjunc.sorted.bam ${outdir}/${id}.Subjunc.bam && samtools index ${outdir}/${id}.Subjunc.sorted.bam ${outdir}/${id}.Subjunc.sorted.bam.bai "
done >subjunc.sh
# 运行
nohup sh subjunc.sh >subjunc.log &
SAM/BAM格式
image.png
image.png
image.png
-
统计比对结果
image.png
# 进入比对文件夹
cd ~/project/Human-16-Asthma-Trans/Mapping/Hisat
# 单个样本
samtools flagstat -@ 3 SRR1039510.Hisat_aln.sorted.bam
ls *.sorted.bam | while read id
do
echo "samtools flagstat -@ 1 ${id} > ${id/bam/flagstat} "
done >flagstat.sh
# 运行
nohup sh flagstat.sh >flagstat.log &
# 质控
multiqc -o ./ *.flagstat
网友评论