bwa
# 建索引
bwa index -a bwtsw ~/hg38.fa
# 100bp长的双端测序reads的比对用bwa的mem算法
bwa mem -t 5 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" /public/biosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38 7E5241.L1_1.fastq 7E5241.L1_2.fastq | samtools sort -@ 5 -o 7E5241.bam -
cat config |while read id
do
arr=($id)
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
bwa mem -t 5 -R "@RG\tID:$sample\tSM:$sample\tLB:WGS\tPL:Illumina" $INDEX $fq1 $fq2 | samtools sort -@ 5 -o $sample.bam -
done
bowtie
# 官网下载索引
# 自己构建
bowtie2-build /home/qmcui/database/reference/hg38/hg38.fa hg38 1>hg38.bowtie2.log 2>hg38.bowtie2.error &
# Total time for backward call to driver() for mirror index: 01:18:29
# 增加线程数
bowtie2-build /home/qmcui/database/reference/hg38/hg38.fa hg38 --threads 5
# Total time for backward call to driver() for mirror index: 00:21:55
# 单端比对
$bin_bowtie2 -p 5 -x $bowtie2_index -U $id | samtools sort -O bam -@ 5 -o - > ${sample}.bam
# 双端序列比对
bowtie2 -p 6 -x /public/reference/index/bowtie/hg38 -1 /home/qmcui/project/7.CHIP/2.val.fq/K_K4Me_clean_1_val_1.fq.gz -2 /home/qmcui/project/7.CHIP/2.val.fq/K_K4Me_clean_2_val_2.fq.gz | samtools sort -O bam -@ 6 -o - >K_K4Me.sort.bam
salmon
直接从fastq到达特异性表达矩阵
# wget ftp://ftp.ensembl.org/pub/release-91/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz
# 注意,其他软件都是利用参考基因组的fa文件建立索引,这个软件用的cdna的数据建立索引
# 建索引
./bin/salmon index -t transcripts.fa -i transcripts_index
# 4min结束
# 双端比对
salmon quant -i transcripts_index -l <LIBTYPE> -1 reads1.fq -2 reads2.fq -o transcripts_quant
# 单端比对
salmon quant -i transcripts_index -l <LIBTYPE> -r reads.fq -o transcripts_quant
# 参考:https://salmon.readthedocs.io/en/latest/salmon.html#using-salmon
subread
# 建索引
subread-build index ../../genome/hg38/hg38.fa -o /home/qmcui/database/index/subread/hg38
subjunc -T 2 -i /home/qmcui/database/index/subread/hg38 -r ${id}_1_val_1.fq -R ${id}_2_val_2.fq -o ${id}.subjunc.sam
备注:需要对应建索引的软件版本,版本不同会报错!
hisat2
# hisat官网
http://ccb.jhu.edu/software/hisat2/index.shtml
# 2018.11.30
# 直接官网下载,不用自己构建
wget -c ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg38.tar.gz
# 自己建立
hisat2-build –p 4 hg38.fa hg38
# 比对
# hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | –sra-acc <SRA accession number>} [-S <hit>]
# https://ccb.jhu.edu/software/hisat2/manual.shtml
hisat2 -p 16 -x ./genome_tran -1 SRR534293_1.fastq -2 SRR534293_2.fastq –S SRR534293.sam
HISAT2建立索引时,就应该把转录组信息加进去。
HISAT2提供两个Python脚本将GTF文件转换成hisat2-build能使用的文件:
extract_exons.py Homo_sapiens.GRCh38.83.chr.gtf > genome.exon
extract_splice_sites.py Homo_sapiens.GRCh38.83.chr.gtf > genome.ss
此外,HISAT2还支持将SNP信息加入到索引中,这样比对的时候就可以考虑SNP的情况。这仍然需要将SNP文件转换成hisat2-build能使用的文件:
extract_snps.py snp142Common.txt > genome.snp
最后,将基因组、转录组、SNP建立索引:
hisat2-build -p 4 genome.fa --snp genome.snp --ss genome.ss --exon genome.exon genome_snp_tran
官网下载展示
压缩包有make_grch38_tran.sh文件,详细记录了创建索引的过程。
$ cat make_hg38.sh
#!/bin/sh
#
# Downloads sequence for the HG38 version of H. spiens (human) from
# UCSC.
#
# The base files, named ??.fa.gz
#
# By default, this script builds and index for just the base files,
# since alignments to those sequences are the most useful. To change
# which categories are built by this script, edit the CHRS_TO_INDEX
# variable below.
#
UCSC_HG38_BASE=http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips
F=hg38.chromFa.tar.gz
get() {
file=$1
if ! wget --version >/dev/null 2>/dev/null ; then
if ! curl --version >/dev/null 2>/dev/null ; then
echo "Please install wget or curl somewhere in your PATH"
exit 1
fi
curl -o `basename $1` $1
return $?
else
wget $1
return $?
fi
}
HISAT2_BUILD_EXE=./hisat2-build
if [ ! -x "$HISAT2_BUILD_EXE" ] ; then
if ! which hisat2-build ; then
echo "Could not find hisat2-build in current directory or in PATH"
exit 1
else
HISAT2_BUILD_EXE=`which hisat2-build`
fi
fi
rm -f genome.fa
get ${UCSC_HG38_BASE}/$F || (echo "Error getting $F" && exit 1)
tar xvzfO $F > genome.fa || (echo "Error unzipping $F" && exit 1)
rm $F
CMD="${HISAT2_BUILD_EXE} genome.fa genome"
echo Running $CMD
if $CMD ; then
echo "genome index built; you may remove fasta files"
else
echo "Index building failed; see error message"
fi
salmon
软件版本older,也会报错!跟subread一样。
# 查看可用的命令
# Salmon v0.9.1
salmon -h
# 查看帮助文档之Salmon's quasi-mapping-based mode
salmon --no-version-check quant --help-reads
# 查看帮助文档之Salmon's alignment-based mode
salmon --no-version-check quant --help-alignment
wget ftp://ftp.ensembl.org/pub/release-91/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz
wget ftp://ftp.ensembl.org/pub/release-91/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.abinitio.fa.gz
wget ftp://ftp.ensembl.org/pub/release-91/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz
wget ftp://ftp.ensembl.org/pub/release-91/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.abinitio.fa.gz
salmon index -t /home/jianmingzeng/reference/transcriptome/Homo_sapiens.GRCh38.cdna.all.fa.gz -i hg38_index
salmon index -t /home/jianmingzeng/reference/transcriptome/Mus_musculus.GRCm38.cdna.all.fa.gz -i mm10_index