走进转录组

STAR序列比对(测试支原体污染序列)

2021-04-09  本文已影响0人  像鸟一样飞过你的高山

转录组测序发现唯一比对率特别低,进一步将未比对到参考基因组的序列提出来比对到NCBI NT库。发现多数序列比对到了Mycoplasma hyorhinis,怀疑是支原体污染,所以我们准备下载Mycoplasma hyorhinis序列参考基因组文件重新比对。

STAR比对结果统计
进入NCBI网站,搜索Mycoplasma hyorhinis(https://www.ncbi.nlm.nih.gov/genome/?term=Mycoplasma+hyorhinis),下载gff和fa文件。
下载建索引参考序列文件 gff和fa文件

使用cufflinks里面gffread将gff转为gtf文件。

gffread GCF_000383515.1_ASM38351v1_genomic.gff -T -o \
 GCF_000383515.1_ASM38351v1_genomic_Mycoplasma_hyorhinis.gtf

使用STAR软件构建索引。STAR建索引特别吃内存,能把你服务器内存用完,然后报错。此时你就要根据你的内存设置limitGenomeGenerateRAM参数,此处设置比所需内存高一点点,如果低于所需内存也会报错。另外线程数可以设置高一点,内存消耗会因为线程数变多而增加,不过不用担心,并不会成倍增加,40个线程内存消耗也就增加了10%。

mkdir S14_INDEX
STAR  --runMode genomeGenerate --runThreadN 10 --genomeDir S14_INDEX/ \
   --genomeFastaFiles GCF_000383515.1_ASM38351v1_genomic.fna \
   --sjdbGTFfile GCF_000383515.1_ASM38351v1_genomic_Mycoplasma_hyorhinis.gtf \
   --sjdbOverhang 149 #reads长度-1

开始比对

R1=SampleName_R1_001.fastq.gz
outdir=./StepAlign/
gtffile=./star_index/GCF_000383515.1_ASM38351v1_genomic_Mycoplasma_hyorhinis.gtf
STAR_index=./star_index/S14_INDEX
R2=${R1//R1_001.paired.fastq.gz/R2_001.paired.fastq.gz}
R1File=`basename ${R1}`
prefix=${R1File//.R1_001.paired.fastq.gz/}
STAR \
  --genomeDir ${STAR_index} \
  --sjdbGTFfile ${gtffile} \
  --limitBAMsortRAM 40000000000 \
  --runThreadN 12 \
  --limitIObufferSize 500000000 \
  --outFilterType BySJout \
  --outFilterMismatchNmax 999 \
  --outFilterMismatchNoverLmax 0.04 \
  --outFilterMultimapNmax 20 \
  --outFilterMatchNminOverLread 0.66 \
  --outFilterIntronMotifs None \
  --outSJfilterReads All \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMunmapped Within \
  --outSAMstrandField intronMotif \
  --outSAMattrRGline ID:${prefix} SM:${prefix} PL:Illumina \
  --alignSJoverhangMin 8 \
  --alignSJDBoverhangMin 1 \
  --alignIntronMin 20 \
  --alignIntronMax 1000000 \
  --alignMatesGapMax 1000000 \
  --chimSegmentMin 15 \
  --chimJunctionOverhangMin 15 \
  --chimScoreMin 0 \
  --chimScoreDropMax 20 \
  --chimScoreSeparation 10 \
  --chimScoreJunctionNonGTAG -1 \
  --quantMode TranscriptomeSAM \
  --quantTranscriptomeBan IndelSoftclipSingleend \
  --outReadsUnmapped Fastx \
  --readFilesIn ${R1} ${R2} \
  --readFilesCommand zcat \
  --outFileNamePrefix ${outdir}/${prefix}.
echo "finish to anign "
上一篇下一篇

猜你喜欢

热点阅读